%% welfare analysis
tic
start_period = (age_sims_br == 1);
sample_half = max([maxiter_distr - 10000,round(sim_periods/2)]);

start_period = zeros(size(age_sims_br));
start_period(sample_half+1,:) = 1;

end_period   = zeros(size(start_period));

total_agents = sum(start_period(:));

 

mpc_reg_br  = NaN(total_agents, 1);
mpc_reg_re  = NaN(total_agents, 1);

mpc_simple_br  = NaN(total_agents, 1);
mpc_simple_re  = NaN(total_agents, 1);

 

%%%

sorted_wealth_br = NaN(sim_periods - sample_half +1, nsim);
sorted_wealth_re = NaN(sim_periods - sample_half +1, nsim);
sorted_wealth_tr = NaN(sim_periods - sample_half +1, nsim);

for i = 1:sim_periods - sample_half+1
    sorted_wealth_br(i,:) = sort( a_sims_br( sample_half -1 + i,:));
    sorted_wealth_re(i,:) = sort( a_sims_re( sample_half -1 + i,:));
    sorted_wealth_tr(i,:) = sort( a_sims_tr( sample_half -1 + i,:));
    
end

wealth_dist_br = mean(sorted_wealth_br)';
wealth_dist_re = mean(sorted_wealth_re)';
wealth_dist_tr = mean(sorted_wealth_tr)';

wealth_dist_br_all = nanmean(a_sims_br)';
wealth_dist_re_all = nanmean(a_sims_re)';
wealth_dist_tr_all = nanmean(a_sims_tr)';

xi_br = linspace(0,0.15,2);
xi_br = [xi_br, linspace(0.15,max(wealth_dist_br),50)];
xi_br = unique(xi_br);


figure
h = histogram( wealth_dist_br,xi_br,'EdgeColor','b','Normalization','probability','DisplayStyle','stairs');
hold on
wealth_pdf_br = h.Values./sum(h.Values);
xi_br = h.BinEdges;
xi_br = mean([xi_br(1:end-1)',xi_br(2:end)'],2);

figure
h_re = histogram(wealth_dist_re,h.BinEdges,'EdgeColor','r','Normalization','probability','DisplayStyle','stairs');
wealth_pdf_re = h_re.Values./sum(h_re.Values);
hold off

%%%%% Figure III.(a) 
figure
plot(xi_br , wealth_pdf_br, 'Color', [0, 0.447, 0.741],'LineWidth',2.5)
hold on
plot(xi_br, wealth_pdf_re,'-.', 'Color', [0.85, 0.075, 0.075],'LineWidth',2.5)
hold off
%plot(xi_br, [wealth_pdf_re',wealth_pdf_br'],'LineWidth',2.5)
xlim([min(xi_br)-0.5, max(xi_br)])
xlim([min(xi_br)-0.5, 25])
ylim([0 0.23])
xlabel('Asset holding (a)')
ylabel('Probability density')
legend('Bounded Rationality Model','Full Information Model','FontSize',13)
title('Wealth Distribution')
ax = gca;
ax.FontSize = 12;
print -dpsc wealth_distq0.eps


clear h h_re

 

wealth_grid = linspace( 0, 5*amax, stateGrid_size)';

cdf_br = NaN(stateGrid_size,1);
cdf_re = cdf_br;

cdf_br_all = NaN(stateGrid_size,1);
cdf_re_all = cdf_br;

cdf_br_mat = NaN(sim_periods - sample_half +1, stateGrid_size);
cdf_re_mat = NaN(sim_periods - sample_half +1, stateGrid_size);


for i = 1:length(wealth_grid)
    
    cdf_br(i) = sum( wealth_dist_br <= wealth_grid(i))./nsim;
    cdf_re(i) = sum( wealth_dist_re <= wealth_grid(i))./nsim;
    
    cdf_br_all(i) = sum( wealth_dist_br_all <= wealth_grid(i))./total_agents;
    cdf_re_all(i) = sum( wealth_dist_re_all <= wealth_grid(i))./total_agents;
    
    cdf_br_mat(:,i) = sum( a_sims_br( sample_half:end,:) <= wealth_grid(i),2)./nsim;
    cdf_re_mat(:,i) = sum( a_sims_re( sample_half:end,:) <= wealth_grid(i),2)./nsim;
    
end

gini_br_mat = NaN(sim_periods - sample_half,1);
gini_re_mat = NaN(sim_periods - sample_half,1);
gini_tr_mat = NaN(sim_periods - sample_half,1);


for i = 1:sim_periods - sample_half
    
    gini_br_mat(i) = gini( ones(nsim,1), (a_sims_br( i+sample_half-1,:).*(a_sims_br(i+sample_half-1,:)>=0))');
    gini_re_mat(i) = gini( ones(nsim,1), (a_sims_re( i+sample_half-1,:).*(a_sims_re(i+sample_half-1,:)>=0))');
    gini_tr_mat(i) = gini( ones(nsim,1), (a_sims_tr( i+sample_half-1,:).*(a_sims_tr(i+sample_half-1,:)>=0))');
    
end

cdf_br_mean = mean(cdf_br_mat)';
cdf_re_mean = mean(cdf_re_mat)';

gini_br_mean = mean(gini_br_mat);
gini_re_mean = mean(gini_re_mat);
gini_tr_mean = mean(gini_tr_mat);


%[gini( ones(nsim,1), wealth_dist_br), gini( ones(nsim,1), wealth_dist_re)];
%[gini_br_mean, gini_re_mean];

perc99_br = prctile( wealth_dist_br, 99);
perc99_re = prctile( wealth_dist_re, 99);
perc99_tr = prctile( wealth_dist_tr, 99);

perc95_br = prctile( wealth_dist_br, 95);
perc95_re = prctile( wealth_dist_re, 95);
perc95_tr = prctile( wealth_dist_tr, 95);

perc90_br = prctile( wealth_dist_br, 90);
perc90_re = prctile( wealth_dist_re, 90);
perc90_tr = prctile( wealth_dist_tr, 90);

wealthsh_top1_br = sum(wealth_dist_br( wealth_dist_br >= perc99_br))/sum(wealth_dist_br);
wealthsh_top1_re = sum(wealth_dist_re( wealth_dist_re >= perc99_re))/sum(wealth_dist_re);
wealthsh_top1_tr = sum(wealth_dist_tr( wealth_dist_tr >= perc99_tr))/sum(wealth_dist_tr);

wealthsh_top5_br = sum(wealth_dist_br( wealth_dist_br >= perc95_br))/sum(wealth_dist_br);
wealthsh_top5_re = sum(wealth_dist_re( wealth_dist_re >= perc95_re))/sum(wealth_dist_re);
wealthsh_top5_tr = sum(wealth_dist_tr( wealth_dist_tr >= perc95_tr))/sum(wealth_dist_tr);


wealthsh_top10_br = sum(wealth_dist_br( wealth_dist_br >= perc90_br))/sum(wealth_dist_br);
wealthsh_top10_re = sum(wealth_dist_re( wealth_dist_re >= perc90_re))/sum(wealth_dist_re);
wealthsh_top10_tr = sum(wealth_dist_tr( wealth_dist_tr >= perc90_tr))/sum(wealth_dist_tr);



agentsh_zeroa_br = sum(  wealth_dist_br == 0)/nsim;
agentsh_zeroa_re = sum(  wealth_dist_re == 0)/nsim;
agentsh_zeroa_re = sum(  wealth_dist_tr == 0)/nsim;

handtomouthsh_br = sum( wealth_dist_br < 0.25*median(wealth_dist_br))/nsim;
handtomouthsh_re = sum( wealth_dist_re < 0.25*median(wealth_dist_re))/nsim;

handtomouthsh_br = sum( wealth_dist_br <= w_br/6)/nsim;
handtomouthsh_re = sum( wealth_dist_re <= w_br/6)/nsim;
handtomouthsh_tr = sum( wealth_dist_tr <= w_tr/6)/nsim;

nega_br = sum( wealth_dist_br < 0)/nsim;
nega_re = sum( wealth_dist_re < 0)/nsim;


wealth_quint_br = prctile( wealth_dist_br, 100*[0.2, 0.4, 0.6, 0.8]);
wealth_quint_re = prctile( wealth_dist_re, 100*[0.2, 0.4, 0.6, 0.8]);
wealth_quint_tr = prctile( wealth_dist_tr, 100*[0.2, 0.4, 0.6, 0.8]);

wealth_sh_quint_br = NaN( length(wealth_quint_br) + 1 ,1);
wealth_sh_quint_re = NaN( length(wealth_quint_br) + 1  ,1);
wealth_sh_quint_tr = NaN( length(wealth_quint_tr) + 1  ,1);

wealth_sh_quint_br(1) = sum( wealth_dist_br( wealth_dist_br <= wealth_quint_br(1)))/sum(wealth_dist_br);
wealth_sh_quint_re(1) = sum( wealth_dist_re( wealth_dist_re <= wealth_quint_re(1)))/sum(wealth_dist_re);
wealth_sh_quint_tr(1) = sum( wealth_dist_tr( wealth_dist_tr <= wealth_quint_tr(1)))/sum(wealth_dist_tr);


for i = 2:length(wealth_quint_br)
    
    wealth_sh_quint_br(i) = sum( wealth_dist_br( wealth_dist_br <= wealth_quint_br(i)))/sum(wealth_dist_br) - sum(wealth_sh_quint_br(1:i-1)) ;
    wealth_sh_quint_re(i) = sum( wealth_dist_re( wealth_dist_re <= wealth_quint_re(i)))/sum(wealth_dist_re) - sum(wealth_sh_quint_re(1:i-1));
    wealth_sh_quint_tr(i) = sum( wealth_dist_tr( wealth_dist_tr <= wealth_quint_tr(i)))/sum(wealth_dist_tr) - sum(wealth_sh_quint_tr(1:i-1));
    
end

wealth_sh_quint_br(end) = sum( wealth_dist_br( wealth_dist_br > wealth_quint_br(end)))/sum(wealth_dist_br);
wealth_sh_quint_re(end) = sum( wealth_dist_re( wealth_dist_re > wealth_quint_re(end)))/sum(wealth_dist_re);
wealth_sh_quint_tr(end) = sum( wealth_dist_tr( wealth_dist_tr > wealth_quint_tr(end)))/sum(wealth_dist_tr);

toc
%% Aguiar et. al. (2020) style regressions

tic
htm_sims_br    = a_sims_br(sample_half:end,:) <= w_br*yi(sample_half:end,:)/6;
htm_sims_re    = a_sims_re(sample_half:end,:) <= w_re*yi(sample_half:end,:)/6;
htm_sims_tr    = a_sims_tr(sample_half:end,:) <= w_tr*yi(sample_half:end,:)/6;

apc_sims_br    = chat_sims_br(sample_half:end,:)./( rt_br*a_sims_br(sample_half:end,:) + w_br*yi(sample_half:end,:));
apc_sims_re    = cstar_sims_re(sample_half:end,:)./( rt_re*a_sims_re(sample_half:end,:) + w_re*yi(sample_half:end,:));
apc_sims_tr    = cstar_sims_tr(sample_half:end,:)./( rt_tr*a_sims_tr(sample_half:end,:) + w_tr*yi(sample_half:end,:));

apc_sims_br_panel = apc_sims_br(3:end,:);
apc_sims_re_panel = apc_sims_re(3:end,:);
apc_sims_tr_panel = apc_sims_tr(3:end,:);

apc_sims_br_panel( age_sims_br(sample_half+2:end,:) <=2) = NaN; 
apc_sims_re_panel( age_sims_br(sample_half+2:end,:) <=2) = NaN; 
apc_sims_tr_panel( age_sims_br(sample_half+2:end,:) <=2) = NaN; 

htm_sims_br_panel = htm_sims_br(1:end-2,:);
htm_sims_br_panel = htm_sims_br_panel(:);


htm_sims_re_panel = htm_sims_re(1:end-2,:);
htm_sims_re_panel = htm_sims_re_panel(:);

htm_sims_tr_panel = htm_sims_tr(1:end-2,:);
htm_sims_tr_panel = htm_sims_tr_panel(:);

mpc_sims_br_panel = mpc_br2(sample_half+2:end,:);
mpc_sims_re_panel = mpc_re(sample_half+2:end,:);
mpc_sims_tr_panel = mpc_tr(sample_half+2:end,:);

mpc_sims_br_panel( age_sims_br(sample_half+2:end,:) <=2) = NaN; 
mpc_sims_re_panel( age_sims_br(sample_half+2:end,:) <=2) = NaN; 
mpc_sims_tr_panel( age_sims_br(sample_half+2:end,:) <=2) = NaN; 


%%%% Income and consumption growth regression

ygrowth_sims_br =log( (rt_br*a_sims_br(sample_half+2:end,:) + w_br.*yi(sample_half+2:end,:))./(rt_br*a_sims_br(sample_half:end-2,:) + w_br.*yi(sample_half:end-2,:)))/2;
ygrowth_sims_re =log( (rt_re*a_sims_re(sample_half+2:end,:) + w_re.*yi(sample_half+2:end,:))./(rt_re*a_sims_re(sample_half:end-2,:) + w_re.*yi(sample_half:end-2,:)))/2;
ygrowth_sims_tr =log( (rt_tr*a_sims_tr(sample_half+2:end,:) + w_tr.*yi(sample_half+2:end,:))./(rt_tr*a_sims_tr(sample_half:end-2,:) + w_tr.*yi(sample_half:end-2,:)))/2;
 
cgrowth_sims_br =log(chat_sims_br(sample_half+2:end,:)./(chat_sims_br(sample_half:end-2,:)))/2;
cgrowth_sims_re =log(cstar_sims_re(sample_half+2:end,:)./(cstar_sims_re(sample_half:end-2,:)))/2;
cgrowth_sims_tr =log(cstar_sims_tr(sample_half+2:end,:)./(cstar_sims_tr(sample_half:end-2,:)))/2;

cgrowth_sims_br( age_sims_br(sample_half+2:end,:) <=5) = NaN; 
cgrowth_sims_re( age_sims_br(sample_half+2:end,:) <=5) = NaN; 
cgrowth_sims_tr( age_sims_br(sample_half+2:end,:) <=5) = NaN; 

[yhtm_reg_br, bint, ~,~,~] = regress( ygrowth_sims_br(:), [ones(length(apc_sims_br_panel(:)),1), htm_sims_br_panel(:)]);
[yhtm_reg_re, bint, ~,~,~] = regress( ygrowth_sims_re(:), [ones(length(apc_sims_re_panel(:)),1), htm_sims_re_panel(:)]);
[yhtm_reg_tr, bint, ~,~,~] = regress( ygrowth_sims_tr(:), [ones(length(apc_sims_tr_panel(:)),1), htm_sims_tr_panel(:)]);

[chtm_reg_br, bint, ~,~,~] = regress( cgrowth_sims_br(:), [ones(length(apc_sims_br_panel(:)),1), htm_sims_br_panel(:)]);
[chtm_reg_re, bint, ~,~,~] = regress( cgrowth_sims_re(:), [ones(length(apc_sims_re_panel(:)),1), htm_sims_re_panel(:)]);
[chtm_reg_tr, bint, ~,~,~] = regress( cgrowth_sims_tr(:), [ones(length(apc_sims_tr_panel(:)),1), htm_sims_tr_panel(:)]);

[mpc_reg_br, bint, ~,~,~] = regress( mpc_sims_br_panel(:), [ones(length(apc_sims_br_panel(:)),1), htm_sims_br_panel(:)]);
[mpc_reg_re, bint, ~,~,~] = regress( mpc_sims_re_panel(:), [ones(length(apc_sims_re_panel(:)),1), htm_sims_re_panel(:)]);
[mpc_reg_tr, bint, ~,~,~] = regress( mpc_sims_tr_panel(:), [ones(length(apc_sims_tr_panel(:)),1), htm_sims_tr_panel(:)]);

%%% fixed effects regressions

start_period = (age_sims_br == 1);
total_agents = sum(start_period(:));

life_length  = NaN(total_agents,1);
chtm_reg_indiv_br = NaN(total_agents,1);
chtm_reg_indiv_re = NaN(total_agents,1);
chtm_reg_indiv_tr = NaN(total_agents,1);

yhtm_reg_indiv_br = NaN(total_agents,1);
yhtm_reg_indiv_re = NaN(total_agents,1);
yhtm_reg_indiv_tr = NaN(total_agents,1);

frac_htm_status_br = NaN(size(a_sims_br));
frac_htm_status_re = NaN(size(a_sims_re));
frac_htm_status_tr = NaN(size(a_sims_re));

abar_br_sims = NaN(total_agents, 1);
rhoa_br_sims = NaN(total_agents, 1);
abar_re_sims = NaN(total_agents, 1);
rhoa_re_sims = NaN(total_agents, 1);

iter = 0;

for i = 1:nsim
    
    ind_start = find(start_period(:,i) == 1);
    ind_end   = [ind_start(2:end)-1;sim_periods];
    
    for j = 1:length(ind_start)
        
        if ind_start(j) < sample_half
            continue
        end
        
        iter = iter + 1;
        
        life_length(iter)  = age_sims_br(ind_end(j),i);
        
        ctemp         = chat_sims_br( ind_start(j):ind_end(j), i);
        ytemp         = w_br.*yi( ind_start(j):ind_end(j), i);
        cgrowth_temp  = log(ctemp(3:end)./ctemp(1:end-2))/2;
        ygrowth_temp  = log(ytemp(3:end)./ytemp(1:end-2))/2;
        
        htm_temp      = a_sims_br( ind_start(j):ind_end(j), i) <= w_br.*yi( ind_start(j):ind_end(j), i)/6;
        %htm_temp      = a_sims_br( ind_start(j):ind_end(j), i) <= w_br/6;
        %frac_htm_status_br(ind_start(j):ind_end(j),i) = mean(htm_temp);
        
        if life_length(iter) > 5 && sum(htm_temp(1:end-2))> 0 && sum(~htm_temp(1:end-2)) > 0
            [b_temp,~,~,~,~] = regress( cgrowth_temp, [ones(length(cgrowth_temp),1), htm_temp(1:end-2)]);
            chtm_reg_indiv_br(iter) = b_temp(2);
            
            [b_temp,~,~,~,~] = regress( ygrowth_temp, [ones(length(cgrowth_temp),1), htm_temp(1:end-2)]);
            yhtm_reg_indiv_br(iter) = b_temp(2);
            frac_htm_status_br(ind_start(j):ind_end(j),i) = mean(htm_temp);
            
        end
        
        %if  life_length(iter) > 5
        %    frac_htm_status_br(ind_start(j):ind_end(j),i) = mean(htm_temp);
        %end
        
        ctemp         = cstar_sims_re( ind_start(j):ind_end(j), i);
        ytemp         = w_re.*yi( ind_start(j):ind_end(j), i);
        cgrowth_temp  =log(ctemp(3:end)./ctemp(1:end-2))/2;
        ygrowth_temp  = log(ytemp(3:end)./ytemp(1:end-2))/2;
        htm_temp      = a_sims_re( ind_start(j):ind_end(j), i) <= w_re.*yi( ind_start(j):ind_end(j), i)/6;
        %htm_temp      = a_sims_re( ind_start(j):ind_end(j), i) <= w_re/6;
 
        
        
        if life_length(iter) > 10 && sum(htm_temp(1:end-2))> 0 && sum(~htm_temp(1:end-2)) > 0
            [b_temp,~,~,~,~] = regress( cgrowth_temp, [ones(length(cgrowth_temp),1), htm_temp(1:end-2)]);
            chtm_reg_indiv_re(iter) = b_temp(2);
            
            [b_temp,~,~,~,~] = regress( ygrowth_temp, [ones(length(cgrowth_temp),1), htm_temp(1:end-2)]);
            yhtm_reg_indiv_re(iter) = b_temp(2);
            frac_htm_status_re(ind_start(j):ind_end(j),i) = mean(htm_temp);
            
        end
        
        %if  sum(htm_temp) > 2
        %    frac_htm_status_re(ind_start(j):ind_end(j),i) = mean(htm_temp);
        %end
        
        ctemp         = cstar_sims_tr( ind_start(j):ind_end(j), i);
        ytemp         = w_tr.*yi( ind_start(j):ind_end(j), i);
        cgrowth_temp  =log(ctemp(3:end)./ctemp(1:end-2))/2;
        ygrowth_temp  = log(ytemp(3:end)./ytemp(1:end-2))/2;
        htm_temp      = a_sims_tr( ind_start(j):ind_end(j), i) <= w_tr.*yi( ind_start(j):ind_end(j), i)/6;
        
        
        
        if life_length(iter) > 10 && sum(htm_temp(1:end-2))> 0 && sum(~htm_temp(1:end-2)) > 0
            [b_temp,~,~,~,~] = regress( cgrowth_temp, [ones(length(cgrowth_temp),1), htm_temp(1:end-2)]);
            chtm_reg_indiv_tr(iter) = b_temp(2);
            
            [b_temp,~,~,~,~] = regress( ygrowth_temp, [ones(length(cgrowth_temp),1), htm_temp(1:end-2)]);
            yhtm_reg_indiv_tr(iter) = b_temp(2);
            frac_htm_status_tr(ind_start(j):ind_end(j),i) = mean(htm_temp);
            
        end
        
        if  life_length(iter)> 10 &&  sum(htm_temp) > 2
            frac_htm_status_tr(ind_start(j):ind_end(j),i) = mean(htm_temp);
        end
        
        
         
    end
    
end


frac_htm_status_br    = frac_htm_status_br(sample_half:end,:);
frac_htm_br_panel = frac_htm_status_br(1:end-2,:);


frac_htm_status_re    = frac_htm_status_re(sample_half:end,:);

frac_htm_status_tr    = frac_htm_status_tr(sample_half:end,:);

temp =  mean( a_sims_re(1:end-2,:) <= w_re/6);
%temp(temp < 2/nsim ) = NaN; 
frac_htm_re_panel    = repmat(temp, [size(htm_sims_re(1:end-2,:),1),1]);

temp =  mean( a_sims_tr(1:end-2,:) <= w_tr/6);
%temp(temp < 2/nsim ) = NaN; 
frac_htm_tr_panel    = repmat(temp, [size(htm_sims_tr(1:end-2,:),1),1]);


%temp =  mean( htm_sims_tr(1:end-2,:));
%frac_htm_tr_panel    = repmat(temp, [size(htm_sims_tr(1:end-2,:),1),1]);

[yhtm_reg_br2, ~, ~,~,~] = regress( ygrowth_sims_br(:), [ones(length(apc_sims_br_panel(:)),1), htm_sims_br_panel(:), frac_htm_br_panel(:)]);
[yhtm_reg_re2, ~, ~,~,~] = regress( ygrowth_sims_re(:), [ones(length(apc_sims_re_panel(:)),1), htm_sims_re_panel(:), frac_htm_re_panel(:)]);
[yhtm_reg_tr2, ~, ~,~,~] = regress( ygrowth_sims_tr(:), [ones(length(apc_sims_tr_panel(:)),1), htm_sims_tr_panel(:), frac_htm_tr_panel(:)]);

[chtm_reg_br2, ~, ~,~,~] = regress( cgrowth_sims_br(:), [ones(length(apc_sims_br_panel(:)),1), htm_sims_br_panel(:), frac_htm_br_panel(:)]);
[chtm_reg_re2, int_re, ~,~,~] = regress( cgrowth_sims_re(:), [ones(length(apc_sims_re_panel(:)),1), htm_sims_re_panel(:), frac_htm_re_panel(:)]);
[chtm_reg_tr2, ~, ~,~,~] = regress( cgrowth_sims_tr(:), [ones(length(apc_sims_tr_panel(:)),1), htm_sims_tr_panel(:), frac_htm_tr_panel(:)]);

toc
 
%%
%%% Learning traps stats

abar_br_sims = NaN(total_agents, 1);
rhoa_br_sims = NaN(total_agents, 1);
abar_re_sims = NaN(total_agents, 1);
rhoa_re_sims = NaN(total_agents, 1);

iter = 0;

flag = 0;



for i = 1:nsim
    
    ind_start = find(start_period(:,i) == 1);
    ind_end   = [ind_start(2:end)-1;sim_periods];
    
    for j = 1:length(ind_start)
        
        if ind_start(j) < sample_half
            continue
        end
        
        iter = iter + 1;
        
        
        
        life_length(iter)  = age_sims_br(ind_end(j),i);
        
        %{
        if life_length(iter) > 150
            flag = 1;
            break
        end
        %}
        
        %%%%% abar agent-by-agent
        if life_length(iter) > 20
            atemp = a_sims_br(ind_start(j) + floor(life_length(iter)/2) :ind_end(j),i)*(1+r_br) + w_br*yi(ind_start(j) + floor(life_length(iter)/2) :ind_end(j),i);
            if sum(atemp > sigma_l*4) < floor(life_length(iter)/4)
                abar_br_sims(iter) = mean(atemp);
                rhoa_br_sims(iter) = 0;
            else
                
                [rho_a,~,~,~,~] = regress( atemp(2:end), [ones(length(atemp(2:end)),1), atemp(1:end-1)]);
                abar_br_sims(iter) = rho_a(1)/(1 - rho_a(2));
                rhoa_br_sims(iter) = rho_a(2);
            end
            
            if ~isfinite(abar_br_sims(iter))
                flag = 1;
                break
            end
            
            atemp = a_sims_re(ind_start(j) + floor(life_length(iter)/2) :ind_end(j),i)*(1+r_re) + w_re*yi(ind_start(j) + floor(life_length(iter)/2) :ind_end(j),i);
            if sum(atemp > sigma_l*4) < floor(life_length(iter)/4)
                abar_br_sims(iter) = mean(atemp);
                rhoa_br_sims(iter) = 0;
            else
                [rho_a,~,~,~,~] = regress( atemp(2:end), [ones(length(atemp(2:end)),1), atemp(1:end-1)]);
                abar_re_sims(iter) = rho_a(1)/(1 - rho_a(2));
                rhoa_re_sims(iter) = rho_a(2);
            end
            
        end
        
    end
    
    if flag == 1
        break
    end
    
    
end


toc
%%
tic
shat_cross = NaN(sim_periods, nsim);
cross_ind_sims = NaN(sim_periods, nsim);
stateGrid_cross = linspace(b,50,1001);
rw_pol_cross = (rt_br/(1+rt_br))*stateGrid_cross' + 1/(1+rt_br)*w_br;
%rw_pol_cross = min( [stateGrid_cross', rw_pol_cross],[],2);

rw_grid_cross_ind = find(stateGrid_cross > rw_pol_cross(1), 1,'first');

pol_sims_size = 50;
chat_pol_sims = NaN(length(rw_pol_cross), pol_sims_size, nsim); 
shat_pol_sims = NaN(length(rw_pol_cross), pol_sims_size, nsim); 

ya_movavg  = NaN(sim_periods, nsim); 

num_zerocross_sims           = NaN(sim_periods, nsim); 
num_learntraps_sims          = NaN(sim_periods, nsim); 
downwardslope_learntrap_sims = NaN(sim_periods, nsim); 
chat_slopedown_sims          = NaN(sim_periods, nsim); 

ya = yi*w_br + a_sims_br*(1+rt_br);
ya_end   = ya(end,:);

dchat_dy = (min(chat_prime_sims(sample_half+1:end,:),ya(sample_half+1:end,:)) - chat_sims_br(sample_half:end-1,:))./(ya(sample_half+1:end,:) - ya(sample_half:end-1,:));

tic
for i = 1:nsim
    
    shat_cross_temp = NaN( sim_periods, 1);
    cross_ind_temp = NaN(sim_periods,1);
    
    ya_avg_temp = NaN(sim_periods,1); 
    
    for j  =  sample_half + 1:sim_periods
        
        ind_start = find(start_period(1:j,i),1,'last');
        %eta_temp =  eta_sims(ind_start :j-1,i);
        
        ind_temp = ind_start:j;
        ind_temp = ind_temp( sigma_etasq_sims(ind_start:j,i) < 1e10);
        
        eta_temp = eta_sims( ind_temp,i);
        %eta_temp = cstar_brsim(yi(ind_temp,i)*w_br + a_sims_br(ind_temp,i)*(1+r_br)) + eps_etai(ind_temp,i)*0.1*sigma_etasq_sims(1,1);
        [chat_pol,shat_pol] =  gp_update( stateGrid_cross' , yi(ind_temp,i)*w_br + a_sims_br(ind_temp,i)*(1+r_br), eta_temp, sigma_etasq_sims(ind_temp,i), cstar_brsim,sigma_c,psi, cov_fun);
        
        %subplot(1,2,1); plot(stateGrid_cross, [chat_pol, cstar_brsim(stateGrid_cross)']) 
        %hold on 
        %scatter(yi(ind_temp,i)*w_br + a_sims_br(ind_temp,i)*(1+r_br), eta_temp)
        %subplot(1,2,2); plot(stateGrid_cross, shat_pol)

        
        if j >= sim_periods - pol_sims_size + 1
            chat_pol_sims(:, j - sim_periods + pol_sims_size,i) = chat_pol; 
            shat_pol_sims(:, j - sim_periods + pol_sims_size,i) = shat_pol; 
        end
        
        pol_diff = rw_pol_cross - chat_pol; 
        cross_ind_vec = find(diff(sign(pol_diff)))+1;
        num_zerocross_sims(j,i)  = length(cross_ind_vec);
        cross_ind_vec = max( cross_ind_vec, ones(size(cross_ind_vec))*rw_grid_cross_ind); 
        
        num_learntraps_sims(j,i) = sum(shat_pol(cross_ind_vec) <   kappa/(2*W_cc)); 
        
        [~,ya_ind] = min( abs(yi(j,i)*w_br + a_sims_br(j,i)*(1+r_br) - stateGrid_cross));
          
        if chat_pol(ya_ind) < rw_pol_cross(ya_ind)
            
            temp =  chat_pol(ya_ind:end)  - rw_pol_cross(ya_ind:end);
            cross_ind = find(temp > 0, 1,'first')-1;
            
            %if isempty(cross_ind)
            %    error('blah')
            %end
            
            if ~isempty(cross_ind)
                cross_ind_temp(j) = cross_ind + ya_ind;
                
            else
                cross_ind_temp(j) = length(stateGrid_cross);
            end
        else
            temp =  chat_pol(1:ya_ind)  - rw_pol_cross(1:ya_ind);
            cross_ind = find(temp < 0, 1,'last');
            if ~isempty(cross_ind)
                cross_ind_temp(j) = cross_ind;
            else
                cross_ind_temp(j) = 1;
            end
            
            %if isempty(cross_ind)
            %    error('blah neg')
            %end
            
        end
        
        if cross_ind_temp(j) < rw_grid_cross_ind
            %cross_ind_temp(j) = find( chat_pol < stateGrid_cross',1,'first');
            cross_ind_temp(j) =rw_grid_cross_ind;
        end
        
        shat_cross_temp(j) = shat_pol(cross_ind_temp(j));       
        ya_avg_temp(j) = mean( yi(ind_start:j-1,i)*w_br + a_sims_br(ind_start:j-1,i)*(1+r_br));
        
        downwardslope_learntrap_sims(j,i) = (dchat_dy(j- sample_half,i) < 0)*(shat_pol(cross_ind_temp(j)) < kappa/(2*W_cc));
        chat_slopedown_sims(j,i) =  sum( diff(chat_pol(:)) < 0) > 0; 
        
    end
    
    shat_cross(:,i) = shat_cross_temp;
    cross_ind_sims(:,i) = cross_ind_temp;
    ya_movavg(:,i) = ya_avg_temp; 
    i;
end
toc



%cross_ind_sims( isnan(cross_ind_sims)) = length(stateGrid_cross);
ya_cross = NaN(sim_periods, nsim);
ya_cross(sample_half+1:end,:) = stateGrid_cross(cross_ind_sims(sample_half+1:end,:));


ya_re = yi*w_re + a_sims_re*(1+rt_re); 

ya_dist_abs = abs( ya_cross(sample_half+1:end,:) - ya(sample_half+1:end,:));
ya_dist_abs_wholesample = abs( ya_cross - ya);
ya_dist = ya_cross(sample_half+1:end,:) - ya(sample_half+1:end,:) ;

yma_dist_abs = abs( ya_movavg(sample_half+1:end,:) - ya(sample_half+1:end,:));


toc
%%

tic

unconstr_ind = htm_sims_br(2:end,:) == 0;
mpc_temp =  mpc_br(sample_half+1:end,:);
alpha_temp = alpha_star_sims(sample_half+1:end,:);
age_temp = age_sims_br(sample_half+1:end,:);


cross_alpha0 = shat_cross <  kappa/(2*W_cc);
cross_alpha0 = cross_alpha0(sample_half+1:end,:);

mpc_unconstr_crossa0 = mpc_temp( unconstr_ind == 1 & cross_alpha0  == 0);
mpc_unconstr_crossa1 = mpc_temp( unconstr_ind == 1 & cross_alpha0  == 1);

 

mpc_htm            = NaN(sim_periods - sample_half,2); 
mpc_aquint         = NaN(sim_periods - sample_half, 5); 
mpc_trap_htm       = NaN(sim_periods - sample_half,2); 
mpc_trap_aquint    = NaN(sim_periods - sample_half,5); 
mpc_notrap_htm     = NaN(sim_periods - sample_half,2); 
mpc_notrap_aquint  = NaN(sim_periods - sample_half,2); 



chat_temp = chat_sims_br(sample_half+1:end,:); 
cstar_temp = min(cstar_brsim(ya), ya - b);
cstar_temp = cstar_temp(sample_half+1:end,:); 

crw_temp  = min(rw_pol_interp( ya ), ya - b); 
crw_temp = crw_temp(sample_half+1:end,:); 

chat_cstar_diff = NaN( sim_periods - sample_half,5); 
chat_cstar_diff_trap = NaN( sim_periods - sample_half,5); 
chat_cstar_diff_notrap = NaN( sim_periods - sample_half,5); 

chat_cstar_absdiff = NaN( sim_periods - sample_half,5); 
chat_cstar_absdiff_trap = NaN( sim_periods - sample_half,5); 
chat_cstar_absdiff_notrap = NaN( sim_periods - sample_half,5); 

chat_crw_diff = chat_cstar_diff; 
chat_crw_diff_trap = chat_cstar_diff; 
chat_crw_diff_notrap = chat_cstar_diff; 

chat_std_diff = chat_cstar_diff; 
chat_std_diff_trap = chat_cstar_diff; 
chat_std_diff_notrap = chat_cstar_diff; 

frac_intrap_aquint = NaN(sim_periods - sample_half, 5); 
frac_intrap_htm_aquint = NaN(sim_periods - sample_half, 5); 
frac_intrap_nothtm_aquint = NaN(sim_periods - sample_half, 5); 

ya_dist_aquint = NaN( sim_periods - sample_half,2); 
ya_dist_posfrac_aquint = NaN(sim_periods - sample_half, 2); 
ya_dist_trap_aquint = NaN( sim_periods - sample_half,2); 
ya_dist_posfrac_trap_aquint= NaN(sim_periods - sample_half, 2); 
ya_dist_notrap_aquint = NaN( sim_periods - sample_half,2); 
ya_dist_notrap_posfrac_aquint = NaN(sim_periods - sample_half, 2); 

apc_cstar_temp    = cstar_temp./( rt_br*a_sims_br(sample_half+1:end,:) + w_br*yi(sample_half+1:end,:));

 
%%%%%%%%%%% Dynan et. al savings rates analysis 
cdelta_bigposs_aterc_br = NaN(sim_periods -sample_half, 3); 
cstardelta_bigposs_aterc_br = NaN(sim_periods -sample_half, 3); 
cdelta_bigposs_aterc_re = NaN(sim_periods -sample_half, 3); 
cdelta_bigposs_medaterc_br =  NaN(sim_periods -sample_half, 3); 

cdelta_bigposs_aterc_perc_br = NaN(sim_periods -sample_half, 3); 
cstardelta_bigposs_aterc_perc_br = NaN(sim_periods -sample_half, 3); 
cdelta_bigposs_aterc_perc_re = NaN(sim_periods -sample_half, 3); 


ya_bigposs_aterc_br = NaN(sim_periods -sample_half, 3); 
ya_bigposs_aterc_re = NaN(sim_periods -sample_half, 3); 

chatrwdiff_bigposs_aterc_br = NaN(sim_periods -sample_half, 3); 
cstarrwdiff_bigposs_aterc_br = NaN(sim_periods -sample_half, 3); 
cstarrwdiff_bigposs_aterc_re = NaN(sim_periods -sample_half, 3); 

adelta_bigposs_aterc_br     = NaN(sim_periods -sample_half, 3); 
ydelta_bigposs_aterc_br     = NaN(sim_periods -sample_half, 3); 

astardelta_bigposs_aterc_br     = NaN(sim_periods -sample_half, 3); 
ystardelta_bigposs_aterc_br     = NaN(sim_periods -sample_half, 3); 

adelta_bigposs_aterc_re     = NaN(sim_periods -sample_half, 3); 
ydelta_bigposs_aterc_re     = NaN(sim_periods -sample_half, 3); 

ygrowth_bigposs_aterc_br     = NaN(sim_periods -sample_half, 3); 


incdelta_bigposs_aterc_br = NaN(sim_periods - sample_half,3);
incdelta_bigposs_aterc_re = NaN(sim_periods - sample_half,3);
earndelta_bigposs_aterc_br = NaN(sim_periods - sample_half,3);
earndelta_bigposs_aterc_re = NaN(sim_periods - sample_half,3);

expydelta_aquint_br = NaN( sim_periods - sample_half, 5); 
expydelta_aquint_re = NaN( sim_periods - sample_half, 5); 
expydelta_aquint_br2 = NaN( sim_periods - sample_half, 5); 

constr_poss_br = NaN(sim_periods - sample_half, 3); 

cdelta_hist_br = []; 

scutoff = 1.27; %%% This is the parameter needed to simulate the Ganong and Noel income drop exercise 
                %%% 1.27 generates 40% income drop, 1.16 for  30% incomde drop 

for i = 1:sim_periods - sample_half
    
    htm_ind_temp = htm_sims_br(1+i,:) == 1;
    quint_temp = prctile( a_sims_br(i + sample_half,:), 100*[0.2, 0.4, 0.6, 0.8]);
    quint_temp_re = prctile( a_sims_re(i + sample_half,:), 100*[0.2, 0.4, 0.6, 0.8]);

    trap_ind_temp = cross_alpha0(i,:) == 1; 
    
    frac_intrap_aquint(i,1) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) <=quint_temp(1) ));
    frac_intrap_aquint(i,2) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   ) );
    frac_intrap_aquint(i,3) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   ) );
    frac_intrap_aquint(i,4) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   ) );
    frac_intrap_aquint(i,5) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(4) ) ); 
    
    frac_intrap_htm_aquint(i,1) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) <=quint_temp(1) & htm_ind_temp == 1));
    frac_intrap_htm_aquint(i,2) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)  & htm_ind_temp == 1  ) );
    frac_intrap_htm_aquint(i,3) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)  & htm_ind_temp == 1 ) );
    frac_intrap_htm_aquint(i,4) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)  & htm_ind_temp == 1 ) );
    frac_intrap_htm_aquint(i,5) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(4) & htm_ind_temp == 1)  ); 
    
    frac_intrap_nothtm_aquint(i,1) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) <=quint_temp(1) & htm_ind_temp == 0));
    frac_intrap_nothtm_aquint(i,2) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2) & htm_ind_temp == 0  ) );
    frac_intrap_nothtm_aquint(i,3) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3) & htm_ind_temp == 0  ) );
    frac_intrap_nothtm_aquint(i,4) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4) & htm_ind_temp == 0  ) );
    frac_intrap_nothtm_aquint(i,5) =  mean(trap_ind_temp( a_sims_br(i + sample_half,:) > quint_temp(4) & htm_ind_temp == 0) ); 
    
    mpc_htm(i,:) = [nanmean(mpc_temp(i, htm_ind_temp == 1  )), nanmean(mpc_temp(i,htm_ind_temp == 0  ))]; 
    mpc_aquint(i,1) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) <=quint_temp(1)  ));
    mpc_aquint(i,2) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)    ) );
    mpc_aquint(i,3) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)    ) );
    mpc_aquint(i,4) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)    ) );
    mpc_aquint(i,5) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(4)   ) ); 
    
    
    
    mpc_trap_htm(i,:) = [nanmean(mpc_temp(i, htm_ind_temp == 1 & trap_ind_temp)), nanmean(mpc_temp(i,htm_ind_temp == 0 & trap_ind_temp))]; 
    mpc_trap_aquint(i,1) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) <=quint_temp(1) & trap_ind_temp));
    mpc_trap_aquint(i,2) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & trap_ind_temp ) );
    mpc_trap_aquint(i,3) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & trap_ind_temp ) );
    mpc_trap_aquint(i,4) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & trap_ind_temp ) );
    mpc_trap_aquint(i,5) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(4) & trap_ind_temp ) ); 
    
    mpc_notrap_htm(i,:)    = [nanmean(mpc_temp(i, htm_ind_temp == 1 & ~trap_ind_temp)), nanmean(mpc_temp(i,htm_ind_temp == 0 & ~trap_ind_temp))]; 
    mpc_notrap_aquint(i,1) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) <=quint_temp(1) & ~trap_ind_temp));
    mpc_notrap_aquint(i,2) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & ~trap_ind_temp ) );
    mpc_notrap_aquint(i,3) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & ~trap_ind_temp ) );
    mpc_notrap_aquint(i,4) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & ~trap_ind_temp ) );
    mpc_notrap_aquint(i,5) =  nanmean(mpc_temp(i  , a_sims_br(i + sample_half,:) > quint_temp(4) & ~trap_ind_temp ) ); 
    
   
    chat_cstar_diff(i,1)  = mean( log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) )));
    chat_cstar_diff(i,2)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)    )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2) )));
    chat_cstar_diff(i,3)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)    )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3) )));
    chat_cstar_diff(i,4)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)    )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4) )));
    chat_cstar_diff(i,5)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   ))                                                   - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   )));

    chat_cstar_diff_trap(i,1)  = mean( log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1)  & trap_ind_temp )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) & trap_ind_temp )));
    chat_cstar_diff_trap(i,2)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)  & trap_ind_temp)));
    chat_cstar_diff_trap(i,3)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)  & trap_ind_temp)));
    chat_cstar_diff_trap(i,4)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)  & trap_ind_temp)));
    chat_cstar_diff_trap(i,5)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   & trap_ind_temp ))                                                   - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   & trap_ind_temp )));

    chat_cstar_diff_notrap(i,1)  = mean( log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1)  & ~trap_ind_temp )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) & ~trap_ind_temp )));
    chat_cstar_diff_notrap(i,2)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & ~trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)  & ~trap_ind_temp)));
    chat_cstar_diff_notrap(i,3)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & ~trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)  & ~trap_ind_temp)));
    chat_cstar_diff_notrap(i,4)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & ~trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)  & ~trap_ind_temp)));
    chat_cstar_diff_notrap(i,5)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   & ~trap_ind_temp ))                                                   - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   & ~trap_ind_temp )));

    chat_cstar_absdiff(i,1)  = mean( abs( log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) ))));
    chat_cstar_absdiff(i,2)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)    )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2) ))));
    chat_cstar_absdiff(i,3)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)    )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3) ))));
    chat_cstar_absdiff(i,4)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)    )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4) ))));
    chat_cstar_absdiff(i,5)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   ))                                                   - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   ))));

    chat_cstar_absdiff_trap(i,1)  = mean( abs(log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1)  & trap_ind_temp )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) & trap_ind_temp ))));
    chat_cstar_absdiff_trap(i,2)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)  & trap_ind_temp))));
    chat_cstar_absdiff_trap(i,3)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)  & trap_ind_temp))));
    chat_cstar_absdiff_trap(i,4)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)  & trap_ind_temp))));
    chat_cstar_absdiff_trap(i,5)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   & trap_ind_temp ))                                                   - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   & trap_ind_temp ))));

    chat_cstar_absdiff_notrap(i,1)  = mean( abs(log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1)  & ~trap_ind_temp )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) & ~trap_ind_temp ))));
    chat_cstar_absdiff_notrap(i,2)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & ~trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)  & ~trap_ind_temp))));
    chat_cstar_absdiff_notrap(i,3)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & ~trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)  & ~trap_ind_temp))));
    chat_cstar_absdiff_notrap(i,4)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & ~trap_ind_temp  )) - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)  & ~trap_ind_temp))));
    chat_cstar_absdiff_notrap(i,5)  = mean( abs(log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   & ~trap_ind_temp ))                                                   - log(cstar_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   & ~trap_ind_temp ))));

    chat_crw_diff(i,1)  = mean( log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) )) - log(crw_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) )));
    chat_crw_diff(i,2)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)    )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2) )));
    chat_crw_diff(i,3)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)    )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3) )));
    chat_crw_diff(i,4)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)    )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4) )));
    chat_crw_diff(i,5)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   ))                                                   - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   )));

    chat_crw_diff_trap(i,1)  = mean( log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1)  & trap_ind_temp )) - log(crw_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) & trap_ind_temp )));
    chat_crw_diff_trap(i,2)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & trap_ind_temp  )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)  & trap_ind_temp)));
    chat_crw_diff_trap(i,3)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & trap_ind_temp  )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)  & trap_ind_temp)));
    chat_crw_diff_trap(i,4)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & trap_ind_temp  )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)  & trap_ind_temp)));
    chat_crw_diff_trap(i,5)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   & trap_ind_temp ))                                                   - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   & trap_ind_temp )));

    chat_crw_diff_notrap(i,1)  = mean( log(chat_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1)  & ~trap_ind_temp )) - log(crw_temp(i, a_sims_br(i + sample_half,:) <=quint_temp(1) & ~trap_ind_temp )));
    chat_crw_diff_notrap(i,2)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & ~trap_ind_temp  )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)  & ~trap_ind_temp)));
    chat_crw_diff_notrap(i,3)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & ~trap_ind_temp  )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)  & ~trap_ind_temp)));
    chat_crw_diff_notrap(i,4)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & ~trap_ind_temp  )) - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)  & ~trap_ind_temp)));
    chat_crw_diff_notrap(i,5)  = mean( log(chat_temp(i,  a_sims_br(i + sample_half,:) > quint_temp(4)   & ~trap_ind_temp ))                                                   - log(crw_temp(i, a_sims_br(i + sample_half,:) > quint_temp(4)   & ~trap_ind_temp )));

    
    ya_dist_trap_aquint(i,1) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) <=quint_temp(1) & trap_ind_temp));
    ya_dist_trap_aquint(i,2) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & trap_ind_temp ) );
    ya_dist_trap_aquint(i,3) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & trap_ind_temp ) );
    ya_dist_trap_aquint(i,4) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & trap_ind_temp ) );
    ya_dist_trap_aquint(i,5) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(4) & trap_ind_temp ) ); 
    
    ya_dist_notrap_aquint(i,1) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) <=quint_temp(1) & ~trap_ind_temp));
    ya_dist_notrap_aquint(i,2) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)   & ~trap_ind_temp ) );
    ya_dist_notrap_aquint(i,3) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)   & ~trap_ind_temp ) );
    ya_dist_notrap_aquint(i,4) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)   & ~trap_ind_temp ) );
    ya_dist_notrap_aquint(i,5) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(4) & ~trap_ind_temp ) ); 
    
    ya_dist_aquint(i,1) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) <=quint_temp(1)  ));
    ya_dist_aquint(i,2) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)     ) );
    ya_dist_aquint(i,3) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)     ) );
    ya_dist_aquint(i,4) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)     ) );
    ya_dist_aquint(i,5) =  median(ya_dist(i  , a_sims_br(i + sample_half,:) > quint_temp(4)   ) ); 
    
    
   
    
    expydelta_aquint_br(i,1) = mean(  rt_br*ya(i + sample_half, a_sims_br(i+sample_half,:) <=quint_temp(1)) + w_br - (1+rt_br)*chat_sims_br(i+sample_half,a_sims_br(i+sample_half,:) <=quint_temp(1))); 
    expydelta_aquint_br(i,2) = mean(  rt_br*ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(1) & a_sims_br(i+sample_half,:) <=quint_temp(2)) + w_br  - (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(1) & a_sims_br(i+sample_half,:) <=quint_temp(2))); 
    expydelta_aquint_br(i,3) = mean(  rt_br*ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(2) & a_sims_br(i+sample_half,:) <=quint_temp(3)) + w_br  - (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(2) & a_sims_br(i+sample_half,:) <=quint_temp(3))); 
    expydelta_aquint_br(i,4) = mean(  rt_br*ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(3) & a_sims_br(i+sample_half,:) <=quint_temp(4)) + w_br  - (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(3) & a_sims_br(i+sample_half,:) <=quint_temp(4)));
    expydelta_aquint_br(i,5) = mean(  rt_br*ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(4)  ) + w_br  - (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(4)  ));

    quint_one = a_sims_br(i+sample_half,:) <=quint_temp(1) & age_sims_br(i + sample_half, :) > 1;
    
    expydelta_aquint_br2(i,1) = mean(  (1+rt_br)*rw_pol_interp(ya(i + sample_half, a_sims_br(i+sample_half,:) <=quint_temp(1) & age_sims_br(i + sample_half, :) > 1))                                              -  (1+rt_br)*chat_sims_br(i+sample_half,a_sims_br(i+sample_half,:) <=quint_temp(1) & age_sims_br(i + sample_half, :) > 1)); 
    expydelta_aquint_br2(i,2) = mean(  (1+rt_br)*rw_pol_interp(ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(1) & a_sims_br(i+sample_half,:) <=quint_temp(2) & age_sims_br(i + sample_half, :) > 1)) -  (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(1) & a_sims_br(i+sample_half,:) <=quint_temp(2) & age_sims_br(i + sample_half, :) > 1)); 
    expydelta_aquint_br2(i,3) = mean(  (1+rt_br)*rw_pol_interp(ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(2) & a_sims_br(i+sample_half,:) <=quint_temp(3) & age_sims_br(i + sample_half, :) > 1)) -  (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(2) & a_sims_br(i+sample_half,:) <=quint_temp(3) & age_sims_br(i + sample_half, :) > 1)); 
    expydelta_aquint_br2(i,4) = mean(  (1+rt_br)*rw_pol_interp(ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(3) & a_sims_br(i+sample_half,:) <=quint_temp(4) & age_sims_br(i + sample_half, :) > 1)) -  (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(3) & a_sims_br(i+sample_half,:) <=quint_temp(4) & age_sims_br(i + sample_half, :) > 1));
    expydelta_aquint_br2(i,5) = mean(  (1+rt_br)*rw_pol_interp(ya(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(4) & age_sims_br(i + sample_half, :) > 1 ))                                             - (1+rt_br)*chat_sims_br(i + sample_half, a_sims_br(i+sample_half,:) > quint_temp(4)  & age_sims_br(i + sample_half, :) > 1));

   
    
   
    terc_temp    = prctile(    a_sims_br(i + sample_half-1,:), 100*[0.475, 0.75]);
    terc_temp_re = prctile( a_sims_re(i + sample_half-1,:), 100*[0.475, 0.75]);

    temp_ind = yi(i + sample_half-1,:) > scutoff & age_sims_br(i + sample_half,:) ~= 1; 
    
    cdelta_bigposs_aterc_br(i,1) = mean(chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    cdelta_bigposs_aterc_br(i ,2) = mean(chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    cdelta_bigposs_aterc_br(i ,3) = mean(chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    constr_poss_br(i, 1) = mean((chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - ya(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)))==0); 
    constr_poss_br(i , 2) = mean((chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - ya(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)))==0); 
    constr_poss_br(i , 3) = mean((chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - ya(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)))==0); 

    cdelta_hist_br = [chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)), cdelta_hist_br]; 
    
    cdelta_bigposs_aterc_perc_br(i ,1) = mean((chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1))))./mean(chat_sims_br(i+sample_half-2, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    cdelta_bigposs_aterc_perc_br(i ,2) = mean((chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))))./mean(chat_sims_br(i+sample_half-2, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    cdelta_bigposs_aterc_perc_br(i ,3) = mean((chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))))./mean(chat_sims_br(i+sample_half-2, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    
    cdelta_bigposs_medaterc_br(i ,1) = median(chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    cdelta_bigposs_medaterc_br(i ,2) = median(chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    cdelta_bigposs_medaterc_br(i ,3) = median(chat_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    ya_temp1 =  ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) -  (1+rt_br)*(cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1))) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1) ) )); 
    ya_temp2 =  ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - (1+rt_br)*(cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))));
    ya_temp3 =  ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - (1+rt_br)*(cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))));
    
    
    cstardelta_bigposs_aterc_br(i ,1) = mean(cstar_brsim(ya_temp1) - cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1))));
    cstardelta_bigposs_aterc_br(i ,2) = mean(cstar_brsim(ya_temp2) - cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))));
    cstardelta_bigposs_aterc_br(i ,3) = mean(cstar_brsim(ya_temp3) - cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))));
    
    cstardelta_bigposs_aterc_perc_br(i ,1) = mean((cstar_brsim(ya_temp1) - cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1))))./cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1))));
    cstardelta_bigposs_aterc_perc_br(i ,2) = mean((cstar_brsim(ya_temp2) - cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))))./cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))));
    cstardelta_bigposs_aterc_perc_br(i ,3) = mean((cstar_brsim(ya_temp3) - cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))))./cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))));
 
   
    cdelta_bigposs_aterc_re(i ,1) = mean(cstar_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) - cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)));
    cdelta_bigposs_aterc_re(i ,2) = mean(cstar_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ) - cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ));
    cdelta_bigposs_aterc_re(i ,3) = mean(cstar_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) - cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)));
    
    cdelta_bigposs_aterc_perc_re(i,1) = mean((cstar_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) - cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1))))./mean(cstar_sims_re(i+sample_half-2, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)));
    cdelta_bigposs_aterc_perc_re(i ,2) = mean((cstar_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2)) - cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2))))./mean(cstar_sims_re(i+sample_half-2, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2)));
    cdelta_bigposs_aterc_perc_re(i ,3) = mean((cstar_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) - cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2))))./mean(cstar_sims_re(i+sample_half-2, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)));
 
    
    ya_bigposs_aterc_br(i ,1) = mean(ya(i+sample_half - 1, temp_ind & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    ya_bigposs_aterc_br(i ,2) = mean(ya(i+sample_half - 1, temp_ind & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    ya_bigposs_aterc_br(i ,3) = mean(ya(i+sample_half - 1, temp_ind & a_sims_br(i+sample_half-1,:) > terc_temp(2)));

    %ygrowth_bigposs_aterc_br(i + sample_half, 1)     = NaN(sim_periods -sample_half, 3); 
    
    ygrowth_bigposs_aterc_br(i ,1) = mean(diff(ya(i+sample_half - 2 : i+sample_half-1, temp_ind & a_sims_br(i+sample_half-1,:) < terc_temp(1))));
    ygrowth_bigposs_aterc_br(i ,2) = mean(diff(ya(i+sample_half - 2 : i+sample_half-1 , temp_ind & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))));
    ygrowth_bigposs_aterc_br(i ,3) = mean(diff(ya(i+sample_half - 2 : i+sample_half-1, temp_ind & a_sims_br(i+sample_half-1,:) > terc_temp(2))));

    
    chatrwdiff_bigposs_aterc_br(i ,1) = mean(chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1) ) )); 
    chatrwdiff_bigposs_aterc_br(i ,2) = mean(chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))));
    chatrwdiff_bigposs_aterc_br(i ,3) = mean(chat_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))));
    
    cstarrwdiff_bigposs_aterc_br(i ,1) = mean(cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1))) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1) ) )); 
    cstarrwdiff_bigposs_aterc_br(i ,2) = mean(cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))));
    cstarrwdiff_bigposs_aterc_br(i ,3) = mean(cstar_brsim(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))) - rw_pol_interp(ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))));
    
    cstarrwdiff_bigposs_aterc_re(i ,1) = mean(cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) - rw_pol_interp_re(ya_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1))));
    cstarrwdiff_bigposs_aterc_re(i ,2) = mean(cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ) - rw_pol_interp_re(ya_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) )));
    cstarrwdiff_bigposs_aterc_re(i ,3) = mean(cstar_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) - rw_pol_interp_re(ya_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2))));
    
    adelta_bigposs_aterc_br(i ,1) = mean(a_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    adelta_bigposs_aterc_br(i ,2) = mean(a_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    adelta_bigposs_aterc_br(i ,3) = mean(a_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    astardelta_bigposs_aterc_br(i ,1) = mean(astar_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - astar_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    astardelta_bigposs_aterc_br(i ,2) = mean(astar_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - astar_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    astardelta_bigposs_aterc_br(i ,3) = mean(astar_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - astar_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    adelta_bigposs_aterc_re(i ,1) = mean(a_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) - a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)));
    adelta_bigposs_aterc_re(i,2) = mean(a_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ) - a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ));
    adelta_bigposs_aterc_re(i ,3) = mean(a_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) - a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)));
    
    ydelta_bigposs_aterc_br(i ,1) = mean(ya(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    ydelta_bigposs_aterc_br(i ,2) = mean(ya(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    ydelta_bigposs_aterc_br(i ,3) = mean(ya(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - ya(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    ydelta_bigposs_aterc_re(i ,1) = mean(ya_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) - ya_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)));
    ydelta_bigposs_aterc_re(i ,2) = mean(ya_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ) - ya_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ));
    ydelta_bigposs_aterc_re(i ,3) = mean(ya_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) - ya_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)));
    
    earndelta_bigposs_aterc_br(i ,1) = mean(w_br*yi(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    earndelta_bigposs_aterc_br(i,2) = mean(w_br*yi(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    earndelta_bigposs_aterc_br(i ,3) = mean(w_br*yi(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    earndelta_bigposs_aterc_re(i ,1) = mean(w_re*yi(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) - w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)));
    earndelta_bigposs_aterc_re(i ,2) = mean(w_re*yi(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ) - w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ));
    earndelta_bigposs_aterc_re(i ,3) = mean(w_re*yi(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) - w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)));
    
    incdelta_bigposs_aterc_br(i ,1) = mean(w_br*yi(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) + rt_br*a_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) - rt_br*a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)))/mean(w_br*yi(i+sample_half-2, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)) + rt_br*a_sims_br(i+sample_half-2, temp_ind == 1 & a_sims_br(i+sample_half-1,:) < terc_temp(1)));
    incdelta_bigposs_aterc_br(i ,2) = mean(w_br*yi(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) + rt_br*a_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2))- w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) - rt_br*a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)))/mean(w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)) + rt_br*a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(1) & a_sims_br(i+sample_half-1,:) < terc_temp(2)));
    incdelta_bigposs_aterc_br(i ,3) = mean(w_br*yi(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) + rt_br*a_sims_br(i+sample_half, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) - w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2))- rt_br*a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)))/mean(w_br*yi(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)) + rt_br*a_sims_br(i+sample_half-1, temp_ind == 1 & a_sims_br(i+sample_half-1,:) > terc_temp(2)));
    
    incdelta_bigposs_aterc_re(i ,1) = mean(w_re*yi(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) + rt_re*a_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1))- w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) - rt_re*a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)))/mean(w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)) + rt_re*a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) < terc_temp_re(1)));
    incdelta_bigposs_aterc_re(i ,2) = mean(w_re*yi(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ) + rt_re*a_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) )- w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) )-  rt_re*a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ))/mean(w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2) ) +  rt_re*a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(1) & a_sims_re(i+sample_half-1,:) < terc_temp_re(2)));
    incdelta_bigposs_aterc_re(i ,3) = mean(w_re*yi(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) + rt_re*a_sims_re(i+sample_half, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2))  - w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) - rt_re*a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)))/mean(w_re*yi(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)) + rt_re*a_sims_re(i+sample_half-1, temp_ind == 1 & a_sims_re(i+sample_half-1,:) > terc_temp_re(2)));

end



htm_sims_br_fullsample    = a_sims_br  <= w_br.*yi /6;



iter = 0; 

trap_freq_lifetime = NaN(total_agents,1);
trap_first_period  = NaN(total_agents,3); 
trap_Tlifelength  = NaN(total_agents,3); 
trap_prob_outoftrap  = NaN(total_agents,3); 
trap_prob_outoftrap_ever  = NaN(total_agents,3); 
rho_ya            = NaN(total_agents,1); 
rho_ya_trap       = NaN(total_agents,1); 
rho_ya_notrap     = NaN(total_agents,1);

rho_yadist            = NaN(total_agents,1); 
rho_yadist_trap       = NaN(total_agents,1); 
rho_yadist_notrap     = NaN(total_agents,1);

Ucstar        = NaN(total_agents,1); 
Uchat         = NaN(total_agents,1); 
Ucstar_notrap = NaN(total_agents,1); 
Uchat_notrap  = NaN(total_agents,1); 
Ucstar_trap = NaN(total_agents,1); 
Uchat_trap  = NaN(total_agents,1); 

cequiv_star        = Ucstar; 
cequiv_star_notrap = Ucstar; 
cequiv_star_trap = Ucstar; 

cequiv_br        = Ucstar; 
cequiv_br_notrap = Ucstar; 
cequiv_br_trap = Ucstar; 

wealth_quint_init = NaN(total_agents,1); 
wealth_quint_trap = NaN(total_agents,1); 

a_trap  = NaN(total_agents,1);
a_init  = NaN(total_agents,1);

life_length = NaN(total_agents,1); 

lifetime_earn_br  = NaN(total_agents,1); 
lifetime_income_br = NaN(total_agents,1); 
lifetime_earn_re  = NaN(total_agents,1); 
lifetime_income_re = NaN(total_agents,1); 
aend_br = NaN(total_agents,1); 
aend_re = NaN(total_agents,1); 

for i = 1:nsim
    
    ind_start = find(start_period(:,i) == 1);
    ind_end   = [ind_start(2:end)-1;sim_periods];
    
    for j = 1:length(ind_start)
        
        if ind_start(j) < sample_half
            continue
        end
        
        iter = iter + 1;
        
        
        
        life_length_temp   = age_sims_br(ind_end(j),i);
        life_length(iter)  = life_length_temp; 
        
        quint_temp = prctile( a_sims_br(ind_start(j),:), 100*[0.2, 0.4, 0.6, 0.8]);
        if ~isempty( find( a_sims_br(ind_start(j),i) >= quint_temp, 1, 'last') )
            wealth_quint_init(iter) = find( a_sims_br(ind_start(j),i) >= quint_temp, 1, 'last') + 1; 
        else
            wealth_quint_init(iter)  = 1;
        end
        
        a_init(iter) = a_sims_br(ind_start(j),i); 
        
        trap_freq_lifetime(iter) = mean(shat_cross(ind_start(j):ind_end(j),i) <  kappa/(2*W_cc));
        cross_temp = shat_cross(ind_start(j):ind_end(j),i) <  kappa/(2*W_cc);
        htm_ind_temp = htm_sims_br_fullsample(ind_start(j):ind_end(j),i); 
        ind_temp = find(shat_cross(ind_start(j):ind_end(j),i) <  kappa/(2*W_cc),1, 'first');
        
        if ~isempty(ind_temp)
            trap_first_period(iter,1) = ind_temp;
            trap_Tlifelength(iter,1)  = ind_temp/life_length_temp;
            if ind_temp <length(cross_temp)
                trap_prob_outoftrap(iter,1)  = mean( cross_temp(ind_temp+1:end) == 0);
                trap_prob_outoftrap_ever(iter,1)  = sum( cross_temp(ind_temp+1:end) == 0)>0;
                
            end
            
            if htm_ind_temp(ind_temp) == 1
                trap_first_period(iter,2) = ind_temp;
                trap_Tlifelength(iter,2)  = ind_temp/life_length_temp;
                if ind_temp <length(cross_temp)
                    trap_prob_outoftrap(iter,2)  = mean( cross_temp(ind_temp+1:end) == 0);
                    trap_prob_outoftrap_ever(iter,2)  = sum( cross_temp(ind_temp+1:end) == 0)>0;
                    
                end
            else
                trap_first_period(iter,3) = ind_temp;
                trap_Tlifelength(iter,3)  = ind_temp/life_length_temp;
                if ind_temp <length(cross_temp)
                    trap_prob_outoftrap(iter,3)  = mean( cross_temp(ind_temp+1:end) == 0);
                    trap_prob_outoftrap_ever(iter,3)  = sum( cross_temp(ind_temp+1:end) == 0)>0;
                    
                end
            end
            
        end
        
        if ~isempty(ind_temp)

            ctemp_star = NaN(life_length_temp - ind_temp+1,1); 
            atemp_star = ctemp_star; 
            astar0     = a_sims_br(ind_start(j):ind_end(j),i); 
            astar0     = astar0(ind_temp); 
            yi_temp    = yi(ind_start(j):ind_end(j),i); 
            yi_temp    = yi_temp(ind_temp:end); 
            
    
            ctemp_star(1)    = min(cstar_brsim(astar0*(1+rt_br) + w_br*yi_temp(1)), astar0*(1+rt_br) + w_br*yi_temp(1) - b);
            atemp_star(1)    = astar0; 

            for i2 = 2:length(ctemp_star)

                atemp_star(i2)    = w_br*yi_temp(i2-1) + (1+rt_br)*atemp_star(i2-1) - ctemp_star(i2-1);
                ctemp_star(i2)    = min(cstar_brsim(atemp_star(i2)*(1+rt_br) + w_br*yi_temp(i2)), atemp_star(i2)*(1+rt_br) + w_br*yi_temp(i2) - b);
            end
            
            a0_temp =  a_sims_br(ind_start(j):ind_end(j),:); 
            quint_temp = prctile( a0_temp(ind_temp,:), 100*[0.2, 0.4, 0.6, 0.8]);
            a_trap(iter) = a0_temp(ind_temp,i);
            
            if ~isempty( find( a0_temp(ind_temp,i) >= quint_temp, 1, 'last') )
                wealth_quint_trap(iter) = find( a0_temp(ind_temp,i) >= quint_temp, 1, 'last') + 1;
            else
                wealth_quint_trap(iter)  = 1;
            end
            
        end
        
        
        if life_length_temp > 10
            b_temp = regress( ya(ind_start(j)+1:ind_end(j),i), [ones(length(ya(ind_start(j):ind_end(j),i))-1,1),ya(ind_start(j):ind_end(j)-1,i)]);
            rho_ya(iter) = b_temp(2);
            
            b_temp = regress( abs( ya_cross(ind_start(j)+1:ind_end(j),i) - ya(ind_start(j)+1:ind_end(j),i)), [ones(length(ya(ind_start(j):ind_end(j),i))-1,1),abs( ya_cross(ind_start(j):ind_end(j)-1,i) - ya(ind_start(j):ind_end(j)-1,i))]);
            rho_yadist(iter) = b_temp(2);
            
            ya_indiv     =  ya(ind_start(j):ind_end(j),i);
            yadist_indiv =  abs(ya_cross(ind_start(j):ind_end(j),i) - ya(ind_start(j):ind_end(j),i));
            
            if length(ya_indiv) - ind_temp > 10
                b_temp = regress( ya_indiv(ind_temp+1:end), [ones(length(ya_indiv(ind_temp+1:end)),1),ya_indiv(ind_temp:end-1)]);
                rho_ya_trap(iter) = b_temp(2);
                
                b_temp = regress( yadist_indiv(ind_temp+1:end), [ones(length(yadist_indiv(ind_temp+1:end)),1),yadist_indiv(ind_temp:end-1)]);
                rho_yadist_trap(iter) = b_temp(2);
            end
            
            if  ind_temp > 10
                b_temp = regress( ya_indiv(2:ind_temp), [ones(length(ya_indiv(2:ind_temp)),1),ya_indiv(1:ind_temp-1)]);
                rho_ya_notrap(iter) = b_temp(2);
                
                b_temp = regress( yadist_indiv(2:ind_temp), [ones(length(yadist_indiv(2:ind_temp)),1),yadist_indiv(1:ind_temp-1)]);
                rho_yadist_notrap(iter) = b_temp(2);
            end
            
        end
        
        %%%%%%%%%%%% Welfare
        
        ctemp = cstar_sims_br( ind_start(j):ind_end(j), i);
        utemp = util(ctemp, util_type, gam);
        
        Ucstar(iter) = sum(utemp.*(bsxfun(@power, btat, 0:length(utemp)-1)'));
        
        ctemp = cstar_sims_br( ind_start(j):ind_end(j), i);
        utemp = util(ctemp(1:ind_temp-1), util_type, gam);
        
        Ucstar_notrap(iter) = sum(utemp.*(bsxfun(@power, btat, 0:length(utemp)-1)'));
        
        ctemp = chat_sims_br( ind_start(j):ind_end(j), i);
        utemp = util(ctemp, util_type, gam);
        
        Uchat(iter) = sum(utemp.*(bsxfun(@power, btat, 0:length(utemp)-1)'));
        
        ctemp = chat_sims_br( ind_start(j):ind_end(j), i);
        utemp = util(ctemp(1:ind_temp-1), util_type, gam);
        
        Uchat_notrap(iter) = sum(utemp.*(bsxfun(@power, btat, 0:length(utemp)-1)'));
        
        if ~isempty(ind_temp)

            ctemp = chat_sims_br( ind_start(j):ind_end(j), i);
            utemp = util(ctemp(ind_temp:end), util_type, gam);

            Uchat_trap(iter) = sum(utemp.*(bsxfun(@power, btat, 0:length(utemp)-1)'));

            ctemp = ctemp_star;
            utemp = util(ctemp , util_type, gam);

            Ucstar_trap(iter) = sum(utemp.*(bsxfun(@power, btat, 0:length(utemp)-1)'));
        end
        
        bta_mult = (1 - btat.^(life_length_temp+1))/(1-btat);
        if ~isempty(ind_temp)
            bta_mult_notrap = (1 - btat.^(ind_temp))/(1-btat);
            bta_mult_trap   = (1 - btat.^(life_length_temp - ind_temp + 1))/(1-btat);
        else 
            bta_mult_notrap = NaN; 
            bta_mult_trap = NaN; 

        end
        %bta_mult = 1/(1-bta);
        
        if gam == 1
            cequiv_star(iter)  = exp(Ucstar(iter)/bta_mult);
            cequiv_star_notrap(iter) = exp(Ucstar_notrap(iter)/bta_mult_notrap);
            cequiv_star_trap(iter) = exp(Ucstar_trap(iter)/bta_mult_trap);

            cequiv_br(iter)  = exp(Uchat(iter)/bta_mult);
            cequiv_br_notrap(iter) = exp(Uchat_notrap(iter)/bta_mult_notrap);
            cequiv_br_trap(iter) = exp(Uchat_trap(iter)/bta_mult_trap);

        else
            cequiv_star(iter)  = (Ucstar(iter)./bta_mult*(1-gam) + 1).^(1/(1-gam));
            cequiv_star_notrap(iter)  = (Ucstar_notrap(iter)./bta_mult_notrap*(1-gam) + 1).^(1/(1-gam));
            cequiv_star_trap(iter)  = (Ucstar_trap(iter)./bta_mult_trap*(1-gam) + 1).^(1/(1-gam));

            cequiv_br(iter)  = (Uchat(iter)./bta_mult*(1-gam) + 1).^(1/(1-gam));
            cequiv_br_notrap(iter)  = (Uchat_notrap(iter)./bta_mult_notrap*(1-gam) + 1).^(1/(1-gam));
            cequiv_br_trap(iter)  = (Uchat_trap(iter)./bta_mult_trap*(1-gam) + 1).^(1/(1-gam));

        end
        
        %%%%%%%% Lifetime earnings and income
        
%         if iter == 2166 
%             error('blah')
%         end
        disc_temp = (1+rt_br).^((ind_start(j):ind_end(j)) - ind_start(j));
        %disc_temp =1; 
        
        %lifetime_earn_br(j,i) = mean(w_br*yi(ind_start:j,i).*(disc_temp(end:-1:1)'));
        %lifetime_income_br(j,i) = mean((rt_br*a_sims_br(ind_start:j,i) + w_br*yi(ind_start:j,i)).*(disc_temp(end:-1:1)'));
        lifetime_earn_br(iter) = sum(w_br*yi( ind_start(j):ind_end(j), i)./(disc_temp'));
        %lifetime_earn_sims_br(j,i) =   sum(w_br*(yi(ind_start:j,i)-1) );
        lifetime_income_br(iter) = sum((rt_br*a_sims_br(ind_start(j):ind_end(j),i) + w_br*yi(ind_start(j):ind_end(j),i))./(disc_temp'));
        %lifetime_income_br(iter) = a_sims_br(ind_start(j),i) + sum(( w_br*yi( ind_start(j):ind_end(j), i))./(disc_temp'));
        %lifetime_income_br(iter) = a_sims_br(ind_start(j), i).*disc_temp(end) + sum(( w_br*yi( ind_start(j):ind_end(j), i)).*(disc_temp(end:-1:1)'));
        aend_br(iter) = a_sims_br(ind_end(j),i); 
        
        disc_temp = (1+rt_re).^((ind_start(j):ind_end(j)) - ind_start(j));
        %disc_temp =1; 

        %lifetime_earn_re(j,i) = mean(w_re*yi(ind_start:j,i).*(disc_temp(end:-1:1)'));
        %lifetime_income_re(j,i) = mean((rt_re*a_sims_re(ind_start:j,i) + w_re*yi(ind_start:j,i)).*(disc_temp(end:-1:1)'));
        lifetime_earn_re(iter) = sum(w_re*yi( ind_start(j):ind_end(j), i)./(disc_temp'));
        %lifetime_earn_sims_br(j,i) =   sum(w_re*(yi(ind_start:j,i)-1) );
        lifetime_income_re(iter) = sum((rt_re*a_sims_re(ind_start(j):ind_end(j),i) + w_re*yi(ind_start(j):ind_end(j),i))./(disc_temp'));
        %lifetime_income_re(iter) = a_sims_re(ind_start(j),i) + sum(( w_re*yi( ind_start(j):ind_end(j), i))./(disc_temp'));
        %lifetime_income_re(iter) = a_sims_re(ind_start(j),i).*disc_temp(end) + sum( (w_re*yi( ind_start(j):ind_end(j), i)).*(disc_temp(end:-1:1)'));
        aend_re(iter) = a_sims_re(ind_end(j),i); 

        
        
    end
    
    
end


quint_temp = prctile( lifetime_earn_br, 100*linspace(0.1,0.9,9)); 
gini_wealth_lifearn_br = NaN( length(quint_temp) + 1,1); 

ind_temp = find( lifetime_earn_br < quint_temp(1)); 
gini_wealth_lifearn_br(1) = gini(ones(length(ind_temp),1), aend_br(ind_temp)); 
    
for i = 2:length(quint_temp) 
    ind_temp = find( lifetime_earn_br > quint_temp(i-1) & lifetime_earn_br <= quint_temp(i)); 
    gini_wealth_lifearn_br(i) = gini(ones(length(ind_temp),1), aend_br(ind_temp)); 

    
end

ind_temp = find( lifetime_earn_br > quint_temp(end)); 
gini_wealth_lifearn_br(end) = gini(ones(length(ind_temp),1), aend_br(ind_temp)); 

quint_temp = prctile( lifetime_income_br, 100*linspace(0.1,0.9,9)); 
gini_wealth_lifinc_br = NaN( length(quint_temp) + 1,1); 

ind_temp = find( lifetime_income_br < quint_temp(1)); 
gini_wealth_lifinc_br(1) = gini(ones(length(ind_temp),1), aend_br(ind_temp)); 
    
for i = 2:length(quint_temp) 
    ind_temp = find( lifetime_income_br > quint_temp(i-1) & lifetime_income_br <= quint_temp(i)); 
    gini_wealth_lifinc_br(i) = gini(ones(length(ind_temp),1), aend_br(ind_temp)); 

    
end

ind_temp = find( lifetime_income_br > quint_temp(end)); 
gini_wealth_lifinc_br(end) = gini(ones(length(ind_temp),1), aend_br(ind_temp)); 


quint_temp = prctile( lifetime_earn_re, 100*linspace(0.1,0.9,9)); 
gini_wealth_lifearn_re = NaN( length(quint_temp) + 1,1); 

ind_temp = find( lifetime_earn_re < quint_temp(1)); 
gini_wealth_lifearn_re(1) = gini(ones(length(ind_temp),1), aend_re(ind_temp)); 

for i = 2:length(quint_temp) 
    ind_temp = find( lifetime_earn_re > quint_temp(i-1) & lifetime_earn_re <= quint_temp(i)); 
    gini_wealth_lifearn_re(i) = gini(ones(length(ind_temp),1), aend_re(ind_temp)); 

    
end

ind_temp = find( lifetime_earn_re > quint_temp(end)); 
gini_wealth_lifearn_re(end) = gini(ones(length(ind_temp),1), aend_re(ind_temp)); 


quint_temp = prctile( lifetime_income_re, 100*linspace(0.1,0.9,9));  
gini_wealth_lifinc_re = NaN( length(quint_temp) + 1 ,1); 

ind_temp = find( lifetime_income_re < quint_temp(1)); 
gini_wealth_lifinc_re(1) = gini(ones(length(ind_temp),1), aend_re(ind_temp)); 
    
for i = 2:length(quint_temp) 
    ind_temp = find( lifetime_income_re > quint_temp(i-1) & lifetime_income_re <= quint_temp(i)); 
    gini_wealth_lifinc_re(i) = gini(ones(length(ind_temp),1), aend_re(ind_temp)); 

    
end

ind_temp = find( lifetime_income_re > quint_temp(end)); 
gini_wealth_lifinc_re(end) = gini(ones(length(ind_temp),1), aend_re(ind_temp)); 



%cequiv_star2  = exp(nanmean(Ucstar))*(1-btat); 
%cequiv_br2  = exp(nanmean(Uchat))*(1-btat); 

chat_temp = chat_sims_br(sample_half+1:end,:); 
cstar_temp = min(cstar_brsim(ya), ya - b);
cstar_temp = cstar_temp(sample_half+1:end,:); 

alpha_prime_temp = alpha_temp(2:end,:); 
cross_alpha0_lag = cross_alpha0(1:end-1,:); 
alpha_lag_temp   = alpha_temp(1:end-1,:); 
cross_alpha0_prime = real(cross_alpha0(2:end,:)); 
cross_alpha0_prime( age_temp(2:end,:) == 1) = NaN; 
age_temp_lag  = age_temp(1:end-1,:); 

ya_dist_shrink = ya_dist_abs(2:end,:)./ya_dist_abs(1:end-1,:); 
age_temp       = age_sims_br(sample_half+1:end,:); 

%%%%% Regressions
ya_dist_lag = ya_dist_abs(1:end-1,:); 
ya_dist_temp = ya_dist_abs(2:end,:); 
ya_dist_temp( age_temp(2:end,:) == 1) = NaN;

[b_yadist, bint,~,~,stats] = regress( ya_dist_temp(:), [ones(length(ya_dist_temp(:)),1), ya_dist_lag(:)]); 
[b_yadist_trap, bint,~,~,stats] = regress( ya_dist_temp(cross_alpha0_lag(:) ==1), [ones(length(ya_dist_temp(cross_alpha0_lag(:) ==1)),1), ya_dist_lag(cross_alpha0_lag(:) ==1)]); 
[b_yadist_notrap, bint,~,~,stats] = regress( ya_dist_temp(cross_alpha0_lag(:) ==0), [ones(length(ya_dist_temp(cross_alpha0_lag(:) ==0)),1), ya_dist_lag(cross_alpha0_lag(:) ==0)]); 

ya_lag = ya(sample_half +  1:end-1,:); 
ya_temp = ya(sample_half + 2:end,:); 
ya_temp( age_temp(2:end,:) == 1) = NaN;

[b_ya, bint,~,~,stats] = regress( ya_temp(:), [ones(length(ya_temp(:)),1), ya_lag(:)]); 
[b_ya_trap, bint,~,~,stats] = regress( ya_temp(cross_alpha0_lag(:) ==1), [ones(length(ya_temp(cross_alpha0_lag(:) ==1)),1), ya_lag(cross_alpha0_lag(:) ==1)]); 
[b_ya_notrap, bint,~,~,stats] = regress( ya_temp(cross_alpha0_lag(:) ==0), [ones(length(ya_temp(cross_alpha0_lag(:) ==0)),1), ya_lag(cross_alpha0_lag(:) ==0)]); 

ya_cross_prime     = ya_cross(sample_half+2:end,:); 
ya_cross_lag       = ya_cross(sample_half+1:end-1,:); 
ya_cross_prime( age_temp(2:end,:) == 1) = NaN; 

corr_alpha_ybardist_trap = corrcoef( alpha_temp(cross_alpha0(:)==1  & age_temp(:) > 1), ya_dist_abs(cross_alpha0(:) == 1 & age_temp(:) > 1),'rows','pairwise'); 
corr_alpha_ybardist_trap = corr_alpha_ybardist_trap(1,2); 
corr_alpha_ybardist_notrap = corrcoef( alpha_temp(cross_alpha0(:)== 0  & age_temp(:) > 1), ya_dist_abs(cross_alpha0(:) == 0 & age_temp(:) > 1),'rows','pairwise'); 
corr_alpha_ybardist_notrap = corr_alpha_ybardist_notrap(1,2); 


corr_alpha_ybardist_trap_posa = corrcoef( alpha_temp(cross_alpha0(:)==1 & alpha_temp(:) >0 & age_temp(:) > 1), ya_dist_abs(cross_alpha0(:) == 1 & alpha_temp(:) >0 & age_temp(:) > 1),'rows','pairwise'); 
corr_alpha_ybardist_trap_posa = corr_alpha_ybardist_trap_posa(1,2); 
corr_alpha_ybardist_notrap_posa = corrcoef( alpha_temp(cross_alpha0(:)== 0 & alpha_temp(:) >0 & age_temp(:) > 1), ya_dist_abs(cross_alpha0(:) == 0 & alpha_temp(:) >0 & age_temp(:) > 1) ,'rows','pairwise'); 
corr_alpha_ybardist_notrap_posa = corr_alpha_ybardist_notrap_posa(1,2); 

corr_alpha_maydist_trap = corrcoef( alpha_temp(cross_alpha0(:)==1  & age_temp(:) > 1), yma_dist_abs(cross_alpha0(:) == 1 & age_temp(:) > 1),'rows','pairwise'); 
corr_alpha_maydist_trap = corr_alpha_maydist_trap(1,2); 
corr_alpha_maydist_notrap = corrcoef( alpha_temp(cross_alpha0(:)== 0  & age_temp(:) > 1), yma_dist_abs(cross_alpha0(:) == 0 & age_temp(:) > 1),'rows','pairwise'); 
corr_alpha_maydist_notrap = corr_alpha_maydist_notrap(1,2); 

htm_ind_temp = htm_sims_br(2:end,:); 
htm_ind_lag  = htm_ind_temp(1:end-1,:); 

corr_alpha_ybardist_trap_htm = corrcoef( alpha_temp(cross_alpha0(:)==1  & age_temp(:) > 1 & htm_ind_temp(:) == 1 ), ya_dist_abs(cross_alpha0(:) == 1 & age_temp(:) > 1 & htm_ind_temp(:) == 1),'rows','pairwise'); 
corr_alpha_ybardist_trap_htm = corr_alpha_ybardist_trap_htm(1,2); 
corr_alpha_ybardist_notrap_htm = corrcoef( alpha_temp(cross_alpha0(:)== 0  & age_temp(:) > 1 & htm_ind_temp(:) == 1 ), ya_dist_abs(cross_alpha0(:) == 0 & age_temp(:) > 1 & htm_ind_temp(:) == 1),'rows','pairwise'); 
corr_alpha_ybardist_notrap_htm = corr_alpha_ybardist_notrap_htm(1,2); 

corr_alpha_ybardist_trap_nothtm = corrcoef( alpha_temp(cross_alpha0(:)==1  & age_temp(:) > 1 & htm_ind_temp(:) == 0 ), ya_dist_abs(cross_alpha0(:) == 1 & age_temp(:) > 1 & htm_ind_temp(:) == 0),'rows','pairwise'); 
corr_alpha_ybardist_trap_nothtm = corr_alpha_ybardist_trap_nothtm(1,2); 
corr_alpha_ybardist_notrap_nothtm = corrcoef( alpha_temp(cross_alpha0(:)== 0  & age_temp(:) > 1 & htm_ind_temp(:) == 0 ), ya_dist_abs(cross_alpha0(:) == 0 & age_temp(:) > 1 & htm_ind_temp(:) == 0),'rows','pairwise'); 
corr_alpha_ybardist_notrap_nothtm = corr_alpha_ybardist_notrap_nothtm(1,2); 


corr_alpha_ybardist_trap_htm_posa = corrcoef( alpha_temp(cross_alpha0(:)==1 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 1 ), ya_dist_abs(cross_alpha0(:) == 1 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 1 ),'rows','pairwise'); 
corr_alpha_ybardist_trap_htm_posa = corr_alpha_ybardist_trap_htm_posa(1,2); 
corr_alpha_ybardist_notrap_htm_posa = corrcoef( alpha_temp(cross_alpha0(:)== 0 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 1 ), ya_dist_abs(cross_alpha0(:) == 0 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 1 ) ,'rows','pairwise'); 
corr_alpha_ybardist_notrap_htm_posa = corr_alpha_ybardist_notrap_htm_posa(1,2); 

corr_alpha_ybardist_trap_nothtm_posa = corrcoef( alpha_temp(cross_alpha0(:)==1 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 0), ya_dist_abs(cross_alpha0(:) == 1 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 0 ),'rows','pairwise'); 
corr_alpha_ybardist_trap_nothtm_posa = corr_alpha_ybardist_trap_nothtm_posa(1,2); 
corr_alpha_ybardist_notrap_nothtm_posa = corrcoef( alpha_temp(cross_alpha0(:)== 0 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 0 ), ya_dist_abs(cross_alpha0(:) == 0 & alpha_temp(:) >0 & age_temp(:) > 1 & htm_ind_temp(:) == 0 ) ,'rows','pairwise'); 
corr_alpha_ybardist_notrap_nothtm_posa = corr_alpha_ybardist_notrap_nothtm_posa(1,2); 

corr_alpha_maydist_trap_htm = corrcoef( alpha_temp(cross_alpha0(:)==1  & age_temp(:) > 1 & htm_ind_temp(:) == 1), yma_dist_abs(cross_alpha0(:) == 1 & age_temp(:) > 1 & htm_ind_temp(:) == 1),'rows','pairwise'); 
corr_alpha_maydist_trap_htm = corr_alpha_maydist_trap_htm(1,2); 
corr_alpha_maydist_notrap_htm = corrcoef( alpha_temp(cross_alpha0(:)== 0  & age_temp(:) > 1 & htm_ind_temp(:) == 1), yma_dist_abs(cross_alpha0(:) == 0 & age_temp(:) > 1 & htm_ind_temp(:) == 1),'rows','pairwise'); 
corr_alpha_maydist_notrap_htm = corr_alpha_maydist_notrap_htm(1,2); 

corr_alpha_maydist_trap_nothtm   = corrcoef( alpha_temp(cross_alpha0(:)==1  & age_temp(:) > 1 & htm_ind_temp(:) == 0), yma_dist_abs(cross_alpha0(:) == 1 & age_temp(:) > 1 & htm_ind_temp(:) == 0),'rows','pairwise'); 
corr_alpha_maydist_trap_nothtm   = corr_alpha_maydist_trap_nothtm(1,2); 
corr_alpha_maydist_notrap_nothtm = corrcoef( alpha_temp(cross_alpha0(:)== 0  & age_temp(:) > 1 & htm_ind_temp(:) == 0), yma_dist_abs(cross_alpha0(:) == 0 & age_temp(:) > 1 & htm_ind_temp(:) == 0),'rows','pairwise'); 
corr_alpha_maydist_notrap_nothtm = corr_alpha_maydist_notrap_nothtm(1,2); 

toc




%% Ergodic policy functions 

tic
%%% Not Normalized, simple average 

chat_pol_avg = nanmean(nanmean( chat_pol_sims, 3),2); 
chat_pol_avg_trap = NaN(size(chat_pol_avg,1), size(chat_pol_sims,2)); 
chat_pol_avg_notrap = chat_pol_avg_trap; 

for j = 1:size(chat_pol_sims,2)
    
    chat_pol_avg_trap(:,j)   = mean(chat_pol_sims(:,j,cross_alpha0( end + j - size(chat_pol_sims,2),:) == 1),3); 
    chat_pol_avg_notrap(:,j) = mean(chat_pol_sims(:,j,cross_alpha0( end+  j - size(chat_pol_sims,2),:) == 0),3); 

end

chat_pol_avg_trap = mean( chat_pol_avg_trap,2); 
chat_pol_avg_notrap = mean( chat_pol_avg_notrap,2); 

chat_pol_aquint_avg = NaN( size(chat_pol_sims,1), size(chat_pol_sims,2),5); 
chat_pol_htm_avg = NaN( size(chat_pol_sims,1), size(chat_pol_sims,2),2); 
chat_pol_aterc_avg = NaN( size(chat_pol_sims,1), size(chat_pol_sims,2),3);  

%shat_pol_aquint_avg = NaN( size(chat_pol_sims,1), size(chat_pol_sims,2),5); 


for j = 1:size(chat_pol_sims,2)
    
    quint_temp = prctile( a_sims_br(j +  sim_periods - size(chat_pol_sims,2),:), 100*[0.2, 0.4, 0.6, 0.8]);
    htm_ind_temp = htm_sims_br( end +  j  - size(chat_pol_sims,2),:) == 1;

    chat_pol_aquint_avg(:,j,1) =  mean(chat_pol_sims(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    chat_pol_aquint_avg(:,j,2) =  mean(chat_pol_sims(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    chat_pol_aquint_avg(:,j,3) =  mean(chat_pol_sims(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    chat_pol_aquint_avg(:,j,4) =  mean(chat_pol_sims(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    chat_pol_aquint_avg(:,j,5) =  mean(chat_pol_sims(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    %shat_pol_aquint_avg(:,j,1) =  mean(shat_pol_sims(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);

    chat_pol_htm_avg(:,j,1) =  mean(chat_pol_sims(:,j, htm_ind_temp ==1 ),3);
    chat_pol_htm_avg(:,j,2) =  mean(chat_pol_sims(:,j, htm_ind_temp ==0 ),3);

    temp_ind    =        yi(j + sim_periods - size(chat_pol_sims,2),:) > scutoff; 
    aterc_temp  = prctile( a_sims_br(j + sim_periods - size(chat_pol_sims,2),:), 100*[0.3333, 0.666]); 
    
    aterc1_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) < aterc_temp(1); 
    aterc2_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) > aterc_temp(1) & a_sims_br(j + sim_periods - size(chat_pol_sims,2) -1,:) < aterc_temp(2); 
    aterc3_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) > aterc_temp(2); 

    chat_pol_aterc_avg(:,j,1) =  mean(chat_pol_sims(:,j, aterc1_ind_temp ==1 ),3);
    chat_pol_aterc_avg(:,j,2) =  mean(chat_pol_sims(:,j, aterc2_ind_temp ==1 ),3);
    chat_pol_aterc_avg(:,j,3) =  mean(chat_pol_sims(:,j, aterc3_ind_temp ==1 ),3);

end

% plot(stateGrid_cross(301:431) , [nanmean(chat_pol_aquint_avg((301:431),:,5),2) , cstar_brsim(stateGrid_cross(301:431) )'])
%plot(stateGrid_cross , [nanmean(chat_pol_aquint_avg(:,:,5),2) , cstar_brsim(stateGrid_cross)'])
%plot(stateGrid_cross , [nanmean(chat_pol_htm_avg(:,:,2),2) , cstar_brsim(stateGrid_cross)'])



%%
%%% Normalized 

chat_pol_norm        = NaN(size(chat_pol_sims)); 

norm_grid_nodes      = 51; 
chat_pol_norm        = NaN( norm_grid_nodes, size(chat_pol_sims,2), size(chat_pol_sims,3)); 
rw_pol_norm          = chat_pol_norm; 
cstar_pol_norm       = chat_pol_norm; 
constr_norm          = chat_pol_norm; 

chat_pol_norm_trap        = NaN( norm_grid_nodes, size(chat_pol_sims,2), size(chat_pol_sims,3)); 
rw_pol_norm_trap          = chat_pol_norm; 
cstar_pol_norm_trap       = chat_pol_norm; 

chat_pol_norm_notrap        = NaN( norm_grid_nodes, size(chat_pol_sims,2), size(chat_pol_sims,3)); 
rw_pol_norm_notrap          = chat_pol_norm; 
cstar_pol_norm_notrap       = chat_pol_norm; 

for j = 1:size(chat_pol_sims,2)
    for i = 1:nsim
        %cross_ind_temp = cross_ind_sims(j + sim_periods - size(chat_pol_sims,2),i);
        [~,cross_ind_temp] = min( abs(  ya(j + sim_periods - size(chat_pol_sims,2),i) - stateGrid_cross));
        
        if  cross_ind_temp > (norm_grid_nodes - 1)/2 && ( size(chat_pol_sims,1) - cross_ind_temp  > (norm_grid_nodes - 1)/2 )
            
            chat_pol_norm(:,j,i)  = chat_pol_sims( cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2 ,j,i); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            rw_pol_norm(:,j,i)    = rw_pol_interp( stateGrid_cross(cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2)); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            cstar_pol_norm(:,j,i) = cstar_brsim( stateGrid_cross(cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2)); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            constr_norm(:,j,i)    = stateGrid_cross(cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2);
            
            if cross_alpha0( end + j  - size(chat_pol_sims,2),i) == 1
                chat_pol_norm_trap(:,j,i)  = chat_pol_sims( cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2 ,j,i);% - rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
                rw_pol_norm_trap(:,j,i)    = rw_pol_interp( stateGrid_cross(cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2));% - rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
                cstar_pol_norm_trap(:,j,i) = cstar_brsim( stateGrid_cross(cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2));% - rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            end
            
            if cross_alpha0( end + j  - size(chat_pol_sims,2),i) == 0
                chat_pol_norm_notrap(:,j,i)  = chat_pol_sims( cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2 ,j,i);% - rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
                rw_pol_norm_notrap(:,j,i)    = rw_pol_interp( stateGrid_cross(cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2));% - rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
                cstar_pol_norm_notrap(:,j,i) = cstar_brsim( stateGrid_cross(cross_ind_temp - (norm_grid_nodes - 1)/2: cross_ind_temp + (norm_grid_nodes - 1)/2));% - rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            end
            
        else
            
            chat_pol_temp = griddedInterpolant(stateGrid_cross, chat_pol_sims(:,j,i),'linear','linear');
            temp_grid = linspace(stateGrid_cross( cross_ind_temp) - ( stateGrid_cross(2) - stateGrid_cross(1))*25,stateGrid_cross( cross_ind_temp) + ( stateGrid_cross(2) - stateGrid_cross(1))*25,51);
            
            chat_pol_norm(:,j,i)  =  chat_pol_temp( temp_grid);
            rw_pol_norm(:,j,i)    =  rw_pol_interp( temp_grid); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            cstar_pol_norm(:,j,i) =  cstar_brsim( temp_grid ); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            constr_norm(:,j,i)    =  temp_grid;

            if cross_alpha0( end  +  j  - size(chat_pol_sims,2),i) == 1
                chat_pol_norm(:,j,i)  =  chat_pol_temp( temp_grid);
                rw_pol_norm(:,j,i)    =  rw_pol_interp( temp_grid); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
                cstar_pol_norm(:,j,i) =  cstar_brsim( temp_grid ); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            end
            
            if cross_alpha0( end + j  - size(chat_pol_sims,2),i) == 0
                chat_pol_norm(:,j,i)  =  chat_pol_temp( temp_grid);
                rw_pol_norm(:,j,i)    =  rw_pol_interp( temp_grid); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
                cstar_pol_norm(:,j,i) =  cstar_brsim( temp_grid ); %- rw_pol_interp(ya_cross(j + sim_periods - size(chat_pol_sims,2),i));
            end
            
            
        end
    end
end

chat_pol_norm_avg = nanmean(nanmean( chat_pol_norm, 3),2); 
rw_pol_norm_avg   = nanmean(nanmean( rw_pol_norm,3),2); 
cstar_pol_norm_avg   = nanmean(nanmean( cstar_pol_norm,3),2); 

stateGrid_norm = stateGrid_cross(1:norm_grid_nodes) - stateGrid_cross( (norm_grid_nodes  + 1)/2);
% 
% plot(stateGrid_norm, [squeeze(chat_pol_norm(:,1,:))])
% plot(stateGrid_norm, [ chat_pol_norm_avg, rw_pol_norm_avg, cstar_pol_norm_avg])

chat_pol_aquint_norm    = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 
rw_pol_aquint_norm      = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 
cstar_pol_aquint_norm   = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 

chat_pol_htm_norm    = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 
rw_pol_htm_norm      = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 
cstar_pol_htm_norm   = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 
constr_htm_norm      = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 

chat_pol_aterc_norm    = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 
rw_pol_aterc_norm      = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 
cstar_pol_aterc_norm   = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 
constr_aterc_norm      = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),2 ); 

chat_pol_aquint_norm_trap    = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 
rw_pol_aquint_norm_trap      = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 
cstar_pol_aquint_norm_trap   = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 

chat_pol_aquint_norm_notrap    = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 
rw_pol_aquint_norm_notrap      = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 
cstar_pol_aquint_norm_notrap   = NaN(size(chat_pol_norm,1), size(chat_pol_norm,2),5 ); 

for j = 1:size(chat_pol_sims,2)
    
    quint_temp = prctile( a_sims_br( j +  sim_periods - size(chat_pol_sims,2),:), 100*[0.2, 0.4, 0.6, 0.8]);
    htm_ind_temp = htm_sims_br( end + j - size(chat_pol_sims,2),:) == 1;
    

    temp_ind =                    yi(j + sim_periods - size(chat_pol_sims,2),:) > scutoff; 
    aterc_temp  = prctile( a_sims_br(j + sim_periods - size(chat_pol_sims,2),:), 100*[0.475,0.75]); 
    
    aterc1_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) < aterc_temp(1); 
    aterc2_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) > aterc_temp(1) & a_sims_br(j + sim_periods - size(chat_pol_sims,2) -1,:) < aterc_temp(2); 
    aterc3_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) > aterc_temp(2); 

    aterc1bigs_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) < aterc_temp(1) & temp_ind == 1; 
    aterc2bigs_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) > aterc_temp(1) & a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) < aterc_temp(2) & temp_ind == 1; 
    aterc3bigs_ind_temp = a_sims_br(j + sim_periods - size(chat_pol_sims,2) ,:) > aterc_temp(2) & temp_ind == 1; 
    
    chat_pol_aquint_norm(:,j,1) =  nanmean(chat_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    chat_pol_aquint_norm(:,j,2) =  nanmean(chat_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    chat_pol_aquint_norm(:,j,3) =  nanmean(chat_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    chat_pol_aquint_norm(:,j,4) =  nanmean(chat_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    chat_pol_aquint_norm(:,j,5) =  nanmean(chat_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    rw_pol_aquint_norm(:,j,1) =  nanmean(rw_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    rw_pol_aquint_norm(:,j,2) =  nanmean(rw_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    rw_pol_aquint_norm(:,j,3) =  nanmean(rw_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    rw_pol_aquint_norm(:,j,4) =  nanmean(rw_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    rw_pol_aquint_norm(:,j,5) =  nanmean(rw_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    cstar_pol_aquint_norm(:,j,1) =  nanmean(cstar_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    cstar_pol_aquint_norm(:,j,2) =  nanmean(cstar_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    cstar_pol_aquint_norm(:,j,3) =  nanmean(cstar_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    cstar_pol_aquint_norm(:,j,4) =  nanmean(cstar_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    cstar_pol_aquint_norm(:,j,5) =  nanmean(cstar_pol_norm(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    chat_pol_aquint_norm_trap(:,j,1) =  nanmean(chat_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    chat_pol_aquint_norm_trap(:,j,2) =  nanmean(chat_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    chat_pol_aquint_norm_trap(:,j,3) =  nanmean(chat_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    chat_pol_aquint_norm_trap(:,j,4) =  nanmean(chat_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    chat_pol_aquint_norm_trap(:,j,5) =  nanmean(chat_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    rw_pol_aquint_norm_trap(:,j,1) =  nanmean(rw_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    rw_pol_aquint_norm_trap(:,j,2) =  nanmean(rw_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    rw_pol_aquint_norm_trap(:,j,3) =  nanmean(rw_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    rw_pol_aquint_norm_trap(:,j,4) =  nanmean(rw_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    rw_pol_aquint_norm_trap(:,j,5) =  nanmean(rw_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    cstar_pol_aquint_norm_trap(:,j,1) =  nanmean(cstar_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    cstar_pol_aquint_norm_trap(:,j,2) =  nanmean(cstar_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    cstar_pol_aquint_norm_trap(:,j,3) =  nanmean(cstar_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    cstar_pol_aquint_norm_trap(:,j,4) =  nanmean(cstar_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    cstar_pol_aquint_norm_trap(:,j,5) =  nanmean(cstar_pol_norm_trap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    chat_pol_aquint_norm_notrap(:,j,1) =  nanmean(chat_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    chat_pol_aquint_norm_notrap(:,j,2) =  nanmean(chat_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    chat_pol_aquint_norm_notrap(:,j,3) =  nanmean(chat_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    chat_pol_aquint_norm_notrap(:,j,4) =  nanmean(chat_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    chat_pol_aquint_norm_notrap(:,j,5) =  nanmean(chat_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    rw_pol_aquint_norm_notrap(:,j,1) =  nanmean(rw_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    rw_pol_aquint_norm_notrap(:,j,2) =  nanmean(rw_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    rw_pol_aquint_norm_notrap(:,j,3) =  nanmean(rw_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    rw_pol_aquint_norm_notrap(:,j,4) =  nanmean(rw_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    rw_pol_aquint_norm_notrap(:,j,5) =  nanmean(rw_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    cstar_pol_aquint_norm_notrap(:,j,1) =  nanmean(cstar_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <=quint_temp(1) ),3);
    cstar_pol_aquint_norm_notrap(:,j,2) =  nanmean(cstar_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(1) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(2)   ) ,3);
    cstar_pol_aquint_norm_notrap(:,j,3) =  nanmean(cstar_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(2) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(3)   ) ,3);
    cstar_pol_aquint_norm_notrap(:,j,4) =  nanmean(cstar_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(3) &  a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) <= quint_temp(4)   ) ,3);
    cstar_pol_aquint_norm_notrap(:,j,5) =  nanmean(cstar_pol_norm_notrap(:,j, a_sims_br(j + sim_periods - size(chat_pol_sims,2),:) > quint_temp(4) ),3 ); 
    
    chat_pol_htm_norm(:,j,1) =  nanmean(chat_pol_norm(:,j, htm_ind_temp == 1 ),3);
    rw_pol_htm_norm(:,j,1) =  nanmean(rw_pol_norm(:,j, htm_ind_temp == 1 ),3);
    cstar_pol_htm_norm(:,j,1) =  nanmean(cstar_pol_norm(:,j, htm_ind_temp == 1 ),3);
    constr_htm_norm(:,j,1)   =   nanmean(constr_norm(:,j, htm_ind_temp == 1),3); 
    
    chat_pol_htm_norm(:,j,2) =  nanmean(chat_pol_norm(:,j, htm_ind_temp == 0 ),3);
    rw_pol_htm_norm(:,j,2) =  nanmean(rw_pol_norm(:,j, htm_ind_temp == 0 ),3);
    cstar_pol_htm_norm(:,j,2) =  nanmean(cstar_pol_norm(:,j, htm_ind_temp == 0 ),3);
    constr_htm_norm(:,j,2)   =   nanmean(constr_norm(:,j, htm_ind_temp == 0),3); 

    chat_pol_aterc_norm(:,j,1)  =  nanmean(chat_pol_norm(:,j, aterc1_ind_temp == 1 ),3);
    rw_pol_aterc_norm(:,j,1)    =  nanmean(rw_pol_norm(:,j, aterc1_ind_temp == 1 ),3);
    cstar_pol_aterc_norm(:,j,1) =  nanmean(cstar_pol_norm(:,j, aterc1_ind_temp == 1 ),3);
    constr_aterc_norm(:,j,1)    =   nanmean(constr_norm(:,j, aterc1_ind_temp == 1),3); 
    
    chat_pol_aterc_norm(:,j,2)  =  nanmean(chat_pol_norm(:,j, aterc2_ind_temp == 1 ),3);
    rw_pol_aterc_norm(:,j,2)    =  nanmean(rw_pol_norm(:,j, aterc2_ind_temp == 1 ),3);
    cstar_pol_aterc_norm(:,j,2) =  nanmean(cstar_pol_norm(:,j, aterc2_ind_temp == 1 ),3);
    constr_aterc_norm(:,j,2)    =  nanmean(constr_norm(:,j, aterc2_ind_temp == 1),3); 
    
    chat_pol_aterc_norm(:,j,3)  =  nanmean(chat_pol_norm(:,j, aterc3_ind_temp == 1 ),3);
    rw_pol_aterc_norm(:,j,3)    =  nanmean(rw_pol_norm(:,j, aterc3_ind_temp == 1 ),3);
    cstar_pol_aterc_norm(:,j,3) =  nanmean(cstar_pol_norm(:,j, aterc3_ind_temp == 1 ),3);
    constr_aterc_norm(:,j,3)    =  nanmean(constr_norm(:,j, aterc3_ind_temp == 1),3); 
end

ya_end_ind = NaN(nsim,1); 

for i = 1:nsim
      
    [~,ya_end_ind(i)] = min( abs( ya_end(i) - stateGrid_cross)); 
    
    
end

chat_polpredict = diag( squeeze(  chat_pol_sims(ya_end_ind ,end,:)) );
diff_temp2 =  cstar_brsim( ya_end)'  - chat_polpredict;

chat_polpredict = min(  [chat_polpredict, ya_end' - b],[],2); 
diff_temp  = chat_sims_br(end,:) - chat_polpredict';

diff_temp3 =  min([cstar_brsim( ya_end)', ya_end' - b],[],2)  - chat_polpredict;

avg_poldiff = NaN( size(chat_pol_aquint_norm,1), 5); 
avg_poldiff(:,1) = nanmean(chat_pol_aquint_norm(:,:,1),2) - nanmean(cstar_pol_aquint_norm(:,:,1),2);
avg_poldiff(:,2) = nanmean(chat_pol_aquint_norm(:,:,2),2) - nanmean(cstar_pol_aquint_norm(:,:,2),2);
avg_poldiff(:,3) = nanmean(chat_pol_aquint_norm(:,:,3),2) - nanmean(cstar_pol_aquint_norm(:,:,3),2);
avg_poldiff(:,4) = nanmean(chat_pol_aquint_norm(:,:,4),2) - nanmean(cstar_pol_aquint_norm(:,:,4),2);
avg_poldiff(:,5) = nanmean(chat_pol_aquint_norm(:,:,5),2) - nanmean(cstar_pol_aquint_norm(:,:,5),2);

clear chat_pol_sims shat_pol_sims

%% Ergodic Policy pictures
 


%%%% Figure 5.(a) 
figure 
plot(stateGrid_norm,nanmean(chat_pol_htm_norm(:,:,1),2), 'LineWidth', 2.5); hold on
plot(stateGrid_norm,nanmean(cstar_pol_htm_norm(:,:,1),2),'--r', 'LineWidth', 2.5);  
plot(stateGrid_norm,nanmean(rw_pol_htm_norm(:,:,1),2), '-.k', 'LineWidth', 1.5);  
plot( stateGrid_norm,nanmean(constr_htm_norm(:,:,1),2),':k','LineWidth',1.5)
ax = gca;
ax.FontSize = 15;
xlabel('Cash-on-hand (Deviation from $y_i$)', 'interpreter','latex')
ylabel('Consumption (c)')
ylim([0.7 1.8])
legend('Bounded rational: $\int \hat c_i(y)$', 'Full info: $c^*(y)$','RW-wealth policy: $c^{RW}(y)$','Borrowing Constraint: $y$','location','Southeast','interpreter','latex','FontSize',15)

print -depsc norm_policy_htm.eps

%%%% Figure 5.(b)
figure 
plot(stateGrid_norm,nanmean(nanmean(chat_pol_aquint_norm(:,:,5),3),2), 'LineWidth', 2.5); hold on
plot(stateGrid_norm,nanmean(nanmean(cstar_pol_aquint_norm(:,:,5),3),2),'--r', 'LineWidth', 2.5);
plot(stateGrid_norm,nanmean(nanmean(rw_pol_aquint_norm(:,:,5),2),3), '-.k', 'LineWidth', 1.5);
ax = gca;
ax.FontSize = 15;
xlabel('Cash-on-hand (Deviation from $y_i$)', 'interpreter','latex')
ylabel('Consumption (c)')
%xlim([0 30])
legend('Bounded rational: $\int \hat c_i(y)$', 'Full info: $c^*(y)$','RW-wealth policy: $c^{RW}(y)$','location','Southeast','interpreter','latex','FontSize',15)

print -depsc norm_policy_q5.eps

%%%% Figure H.1 
figure 
plot(stateGrid_cross, nanmean(chat_pol_aterc_avg(:,:,1),2) + 0.1, '-b','LineWidth',2)
hold on 
plot(stateGrid_cross, cstar_brsim(stateGrid_cross )','-r','LineWidth',2)
plot(stateGrid_cross, rw_pol_cross,'-.k','LineWidth',1.25)
scatter( stateGrid_cross(39),nanmean(chat_pol_aterc_avg(39,:,1),2) + 0.1,200,'b','filled')
scatter( stateGrid_cross(39),cstar_brsim(stateGrid_cross(39)),200,'r','filled')
scatter( stateGrid_cross(29),stateGrid_cross(29),200,'b')
scatter( stateGrid_cross(38),cstar_brsim(stateGrid_cross(38)),200,'r')
plot(stateGrid_cross, stateGrid_cross' - b,':k','LineWidth',1.25)
ax = gca;
ax.FontSize = 10;
xlabel('Cash on hand (y)')
ylabel('Consumption (c)')
xlim([0.5 4])
ylim([0.8 1.75])
legend('Avg policy estimate: $\int \widehat c_i(y) |$ Low $a_{t-1}$ ', 'Full info: $c^*(y)$','RW-y policy: $c^{RW}(y)$','Time-$t$ action: $c_t$','Time-$t$ opt. action: $ c_t^*$','Time-$t+1$ action: $c_{t+1}$','Time-$t$ opt. action: $ c_{t+1}^*$','location','Southeast','interpreter','latex','FontSize',12)

print -depsc avg_policy_aterc1.eps


%% ybar_i disribution 

xi_ybar = linspace(min(ya_cross(:)), max(ya_cross(:)), 45); 
xi_ybar = sort( [xi_ybar, mean(xi_ybar(1:2))]);
%xi_ybar = [linspace(xi_ybar(1), mean(xi_ybar(1:2)),3), xi_ybar(2:end)];

h = histogram( ya_cross(:), xi_ybar,'EdgeColor', 'b','Normalization','probability','DisplayStyle','stairs')
ybari_pdf = h.Values./sum(h.Values); 
xi_ybar = h.BinEdges; 
xi_ybar = mean( [xi_ybar(1:end-1)', xi_ybar(2:end)'],2); 


fi_cross_ind= find( cstar_brsim(stateGrid_cross)' - rw_pol_cross > 0,1,'first');
if isempty(fi_cross_ind)
    fi_cross_ind =1;
end

%%%%% Figure III.(b) 
figure
ax = axes;
plot(xi_ybar, ybari_pdf,'Color',[0, 0.4470, 0.7410],'LineWidth',2.5)
hold on
xline(stateGrid_cross(fi_cross_ind),'-.','color',[0.8500, 0.075, 0.075],'LineWidth',2.5)
xlim([min(xi_ybar)-0.5, 25])
ylim([0 0.35])
xlabel('Cash-on-hand (y)')
ylabel('Probability density')
%yticks([0,0.05,0.1,0.15,0.2])
legend('Bounded Rationality Model','Full Information Model','FontSize',15)

%title('Wealth Distribution')
%ax = gca;
ax.FontSize = 12;

print -depsc ybari_dist.eps
 
%}

%% MPCs
mean_mpc_br_sims = mean(mpc_br,2);
mean_mpc_re_sims = mean(mpc_re,2);
mean_mpc_br_sims2 = mean(mpc_br2,2);
mean_mpc_re_sims2 = mean(mpc_re2,2);
mean_mpc_br_sims3 = mean(mpc_br3,2);
mean_mpc_re_sims3 = mean(mpc_re3,2);
mean_mpc_br_sims4 = mean(mpc_br4,2);
mean_mpc_re_sims4 = mean(mpc_re4,2);

median_mpc_br_sims = median(mpc_br,2);
median_mpc_re_sims = median(mpc_re,2);
median_mpc_br_sims2 = median(mpc_br2,2);
median_mpc_re_sims2 = median(mpc_re2,2);
median_mpc_br_sims3 = median(mpc_br3,2);
median_mpc_re_sims3 = median(mpc_re3,2);
median_mpc_br_sims4 = median(mpc_br4,2);
median_mpc_re_sims4 = median(mpc_re4,2);


mean_mpc_tr_sims = mean(mpc_tr,2);
mean_mpc_tr_sims2 = mean(mpc_tr2,2);
mean_mpc_tr_sims3 = mean(mpc_tr3,2);
mean_mpc_tr_sims4 = mean(mpc_tr4,2);

median_mpc_tr_sims = median(mpc_tr,2);
median_mpc_tr_sims2 = median(mpc_tr2,2);
median_mpc_tr_sims3 = median(mpc_tr3,2);
median_mpc_tr_sims4 = median(mpc_tr4,2);

std_mpc_br_sims = std(mpc_br,1,2);
std_mpc_re_sims = std(mpc_re,1,2);
std_mpc_tr_sims = std(mpc_tr,1,2);


prct5_mpc_br_sims = prctile(mpc_br,5,2);
prct5_mpc_re_sims = prctile(mpc_re,5,2);
prct10_mpc_br_sims = prctile(mpc_br,10,2);
prct10_mpc_re_sims = prctile(mpc_re,10,2);
prct90_mpc_br_sims = prctile(mpc_br,90,2);
prct90_mpc_re_sims = prctile(mpc_re,90,2);
prct95_mpc_br_sims = prctile(mpc_br,95,2);
prct95_mpc_re_sims = prctile(mpc_re,95,2);

prct5_mpc_tr_sims = prctile(mpc_tr,5,2);
prct10_mpc_tr_sims = prctile(mpc_tr,10,2);
prct90_mpc_tr_sims = prctile(mpc_tr,90,2);
prct95_mpc_tr_sims = prctile(mpc_tr,95,2);

mean_mpc_br = mean( mean_mpc_br_sims(sample_half:end));
mean_mpc_re = mean( mean_mpc_re_sims(sample_half:end));
mean_mpc_br2 = mean( mean_mpc_br_sims2(sample_half:end));
mean_mpc_re2 = mean( mean_mpc_re_sims2(sample_half:end));
mean_mpc_br3 = mean( mean_mpc_br_sims3(sample_half:end));
mean_mpc_re3 = mean( mean_mpc_re_sims3(sample_half:end));
mean_mpc_br4 = mean( mean_mpc_br_sims4(sample_half:end));
mean_mpc_re4 = mean( mean_mpc_re_sims4(sample_half:end));

mean_mpc_tr = mean( mean_mpc_tr_sims(sample_half:end));
mean_mpc_tr2 = mean( mean_mpc_tr_sims2(sample_half:end));
mean_mpc_tr3 = mean( mean_mpc_tr_sims3(sample_half:end));
mean_mpc_tr4 = mean( mean_mpc_tr_sims4(sample_half:end));

median_mpc_br = mean( median_mpc_br_sims(sample_half:end));
median_mpc_re = mean( median_mpc_re_sims(sample_half:end));
median_mpc_br2 = mean( median_mpc_br_sims2(sample_half:end));
median_mpc_re2 = mean( median_mpc_re_sims2(sample_half:end));
median_mpc_br3 = mean( median_mpc_br_sims3(sample_half:end));
median_mpc_re3 = mean( median_mpc_re_sims3(sample_half:end));
median_mpc_br4 = mean( median_mpc_br_sims4(sample_half:end));
median_mpc_re4 = mean( median_mpc_re_sims4(sample_half:end));

median_mpc_tr = mean( median_mpc_tr_sims(sample_half:end));
median_mpc_tr2 = mean( median_mpc_tr_sims2(sample_half:end));
median_mpc_tr3 = mean( median_mpc_tr_sims3(sample_half:end));
median_mpc_tr4 = mean( median_mpc_tr_sims4(sample_half:end));

prct5_mpc_br=mean(prct5_mpc_br_sims(sample_half:end));
prct5_mpc_re=mean(prct5_mpc_re_sims(sample_half:end));
prct10_mpc_br=mean(prct10_mpc_br_sims(sample_half:end));
prct10_mpc_re=mean(prct10_mpc_re_sims(sample_half:end));
prct90_mpc_br=mean(prct90_mpc_br_sims(sample_half:end));
prct90_mpc_re=mean(prct90_mpc_re_sims(sample_half:end));
prct95_mpc_br=mean(prct95_mpc_br_sims(sample_half:end));
prct95_mpc_re=mean(prct95_mpc_re_sims(sample_half:end));

prct5_mpc_tr=mean(prct5_mpc_tr_sims(sample_half:end));
prct10_mpc_tr=mean(prct10_mpc_tr_sims(sample_half:end));
prct90_mpc_tr=mean(prct90_mpc_tr_sims(sample_half:end));
prct95_mpc_tr=mean(prct95_mpc_tr_sims(sample_half:end));


std_mpc_br = mean( std_mpc_br_sims(sample_half:end));
std_mpc_re = mean( std_mpc_re_sims(sample_half:end));
std_mpc_tr = mean( std_mpc_tr_sims(sample_half:end));

mean_mpc_br_abovemed_sims = NaN( sim_periods -sample_half,1);
mean_mpc_re_abovemed_sims = NaN( sim_periods -sample_half,1);
mean_mpc_br_nonza_sims = NaN( sim_periods -sample_half,1);
mean_mpc_re_nonza_sims = NaN( sim_periods -sample_half,1);

mpc_br_aquint_sims = NaN(sim_periods - sample_half,5);
mpc_re_aquint_sims = NaN(sim_periods - sample_half,5);

mpc_br_top1_sims  = NaN(sim_periods - sample_half,1);
mpc_br_top10_sims = NaN(sim_periods - sample_half,1);

mpc_re_top1_sims  = NaN(sim_periods - sample_half,1);
mpc_re_top10_sims = NaN(sim_periods - sample_half,1);

mean_mpc_tr_abovemed_sims = NaN( sim_periods -sample_half,1);
mean_mpc_tr_nonza_sims = NaN( sim_periods -sample_half,1);

mpc_tr_aquint_sims = NaN(sim_periods - sample_half,5);

mpc_tr_top1_sims  = NaN(sim_periods - sample_half,1);
mpc_tr_top10_sims = NaN(sim_periods - sample_half,1);


for i = 1:sim_periods - sample_half
    
    mean_mpc_br_abovemed_sims(i) = mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) > median(a_sims_br(i + sample_half,:))));
    mean_mpc_re_abovemed_sims(i) = mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) > median(a_sims_re(i + sample_half,:))));
    mean_mpc_tr_abovemed_sims(i) = mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) > median(a_sims_tr(i + sample_half,:))));
    
    mean_mpc_br_nonza_sims(i) = mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) > 0));
    mean_mpc_re_nonza_sims(i) = mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) > 0));
    mean_mpc_tr_nonza_sims(i) = mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) > 0));
    
    quint_temp = prctile( a_sims_br(i + sample_half,:), 100*[0.2, 0.4, 0.6, 0.8]);
    mpc_br_aquint_sims(i,1) =  mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) <= quint_temp(1)));
    mpc_br_aquint_sims(i,2) =  mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)));
    mpc_br_aquint_sims(i,3) =  mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)));
    mpc_br_aquint_sims(i,4) =  mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)));
    mpc_br_aquint_sims(i,5) =  mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(4)));
    
    quint_temp = prctile( a_sims_re(i + sample_half,:), 100*[0.2, 0.4, 0.6, 0.8]);
    mpc_re_aquint_sims(i,1) =  mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) <= quint_temp(1)));
    mpc_re_aquint_sims(i,2) =  mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) > quint_temp(1) &  a_sims_re(i + sample_half,:) <= quint_temp(2)));
    mpc_re_aquint_sims(i,3) =  mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) > quint_temp(2) &  a_sims_re(i + sample_half,:) <= quint_temp(3)));
    mpc_re_aquint_sims(i,4) =  mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) > quint_temp(3) &  a_sims_re(i + sample_half,:) <= quint_temp(4)));
    mpc_re_aquint_sims(i,5) =  mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) > quint_temp(4)));
    
    quint_temp = prctile( a_sims_tr(i + sample_half,:), 100*[0.2, 0.4, 0.6, 0.8]);
    mpc_tr_aquint_sims(i,1) =  mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) <= quint_temp(1)));
    mpc_tr_aquint_sims(i,2) =  mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) > quint_temp(1) &  a_sims_tr(i + sample_half,:) <= quint_temp(2)));
    mpc_tr_aquint_sims(i,3) =  mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) > quint_temp(2) &  a_sims_tr(i + sample_half,:) <= quint_temp(3)));
    mpc_tr_aquint_sims(i,4) =  mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) > quint_temp(3) &  a_sims_tr(i + sample_half,:) <= quint_temp(4)));
    mpc_tr_aquint_sims(i,5) =  mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) > quint_temp(4)));
    
    prct_temp = prctile(a_sims_br(i + sample_half,:), [90,99]);
    
    mpc_br_top1_sims(i)  =  mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) >= prct_temp(2)));
    mpc_br_top10_sims(i) =  mean(mpc_br(i + sample_half, a_sims_br(i + sample_half,:) >= prct_temp(1)));
    
    prct_temp = prctile(a_sims_re(i + sample_half,:), [90,99]);
    mpc_re_top1_sims(i)  =  mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) >= prct_temp(2)));
    mpc_re_top10_sims(i) =  mean(mpc_re(i + sample_half, a_sims_re(i + sample_half,:) >= prct_temp(1)));
    
    prct_temp = prctile(a_sims_tr(i + sample_half,:), [90,99]);
    mpc_tr_top1_sims(i)  =  mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) >= prct_temp(2)));
    mpc_tr_top10_sims(i) =  mean(mpc_tr(i + sample_half, a_sims_tr(i + sample_half,:) >= prct_temp(1)));
    
end

mean_mpc_br_abovemed = mean(mean_mpc_br_abovemed_sims);
mean_mpc_re_abovemed = mean(mean_mpc_re_abovemed_sims);
mean_mpc_br_nonza = mean(mean_mpc_br_nonza_sims);
mean_mpc_re_nonza = mean(mean_mpc_re_nonza_sims);

mean_mpc_tr_abovemed = mean(mean_mpc_tr_abovemed_sims);
mean_mpc_tr_nonza = mean(mean_mpc_tr_nonza_sims);

 

mpc_htm_br = mpc_br(sample_half:end,:);
mpc_htm_br = mpc_htm_br(htm_sims_br == 1);

mpc_nonhtm_br = mpc_br(sample_half:end,:);
mpc_nonhtm_br = mpc_nonhtm_br(htm_sims_br == 0);
mpc_nonhtm_br2 = mpc_br2(sample_half:end,:);
mpc_nonhtm_br2 = mpc_nonhtm_br2(htm_sims_br == 0);
mpc_nonhtm_br3 = mpc_br3(sample_half:end,:);
mpc_nonhtm_br3 = mpc_nonhtm_br3(htm_sims_br == 0);
mpc_nonhtm_br4 = mpc_br4(sample_half:end,:);
mpc_nonhtm_br4 = mpc_nonhtm_br4(htm_sims_br == 0);

mpc_htm_re = mpc_re(sample_half:end,:);
mpc_htm_re = mpc_htm_re(htm_sims_re == 1);

mpc_nonhtm_re = mpc_re(sample_half:end,:);
mpc_nonhtm_re = mpc_nonhtm_re(htm_sims_re == 0);
mpc_nonhtm_re2 = mpc_re2(sample_half:end,:);
mpc_nonhtm_re2 = mpc_nonhtm_re2(htm_sims_re == 0);
mpc_nonhtm_re3 = mpc_re3(sample_half:end,:);
mpc_nonhtm_re3 = mpc_nonhtm_re3(htm_sims_re == 0);
mpc_nonhtm_re4 = mpc_re4(sample_half:end,:);
mpc_nonhtm_re4 = mpc_nonhtm_re4(htm_sims_re == 0);

mpc_htm_tr = mpc_tr(sample_half:end,:);
mpc_htm_tr = mpc_htm_tr(htm_sims_tr == 1);

mpc_nonhtm_tr = mpc_tr(sample_half:end,:);
mpc_nonhtm_tr = mpc_nonhtm_tr(htm_sims_tr == 0);
mpc_nonhtm_tr2 = mpc_tr2(sample_half:end,:);
mpc_nonhtm_tr2 = mpc_nonhtm_tr2(htm_sims_tr == 0);
mpc_nonhtm_tr3 = mpc_tr3(sample_half:end,:);
mpc_nonhtm_tr3 = mpc_nonhtm_tr3(htm_sims_tr == 0);
mpc_nonhtm_tr4 = mpc_tr4(sample_half:end,:);
mpc_nonhtm_tr4 = mpc_nonhtm_tr4(htm_sims_tr == 0);

%%%% Simple MPC calculations
ya_br_sims = a_sims_br*(1+rt_br) + w_br*yi;
diff_ya_br_sims = diff(ya_br_sims);
diff_ya_br_sims(reset_shocks(2:end,:) < deltax) = NaN;
diff_chat_br_sims = diff(chat_sims_br);

simple_mpc_br_sims  = NaN( size(a_sims_br));
diff_chat_br_sims(age_sims_br(2:end,:) == 1) = NaN; 
simple_mpc_br_sims(2:end,:) = diff_chat_br_sims./diff_ya_br_sims;

median_simple_mpc_br = mean(nanmedian( simple_mpc_br_sims(sample_half:end,:),2));


%%% Reasoning statistics

mean_alpha = mean(mean(alpha_star_sims(sample_half+1:end,:)));
std_alpha  = mean(std(alpha_star_sims(sample_half+1:end,:)));
median_alpha  = mean(median(alpha_star_sims(sample_half+1:end,:)));

alpha_dist = NaN(sim_periods - sample_half, nsim);
alphazero_frac = NaN(sim_periods - sample_half, 1);
avg_alpha_alphaquint = NaN(sim_periods - sample_half, 5);
avg_alpha_aquint = NaN(sim_periods - sample_half, 5);
alphanz_frac_aquint = NaN(sim_periods - sample_half,5);

for i = 1:sim_periods - sample_half
    
    alpha_sort_temp =  sort(alpha_star_sims(sample_half+1,:));
    alpha_dist(i,:) = alpha_sort_temp;
    alphazero_frac(i) = sum(alpha_sort_temp == 0)/nsim;
    
    %quint_temp = prctile( alpha_star_sims(sample_half+1,:), 100*[0.2,0.4,0.6,0.8]);
    avg_alpha_alphaquint(i,1) = mean( alpha_sort_temp(1:nsim/5));
    avg_alpha_alphaquint(i,2) = mean( alpha_sort_temp(nsim/5+1:2*nsim/5));
    avg_alpha_alphaquint(i,3) = mean( alpha_sort_temp(2*nsim/5+1:3*nsim/5));
    avg_alpha_alphaquint(i,4) = mean( alpha_sort_temp(3*nsim/5+1:4*nsim/5));
    avg_alpha_alphaquint(i,5) = mean( alpha_sort_temp(4*nsim/5+1:end));
    
    quint_temp = prctile( a_sims_br(i + sample_half,:), 100*[0.2, 0.4, 0.6, 0.8]);
    avg_alpha_aquint(i,1) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) <= quint_temp(1)));
    avg_alpha_aquint(i,2) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)));
    avg_alpha_aquint(i,3) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)));
    avg_alpha_aquint(i,4) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)));
    avg_alpha_aquint(i,5) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(4)));
    
    alphanz_frac_aquint(i,1) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) <= quint_temp(1)) > 0);
    alphanz_frac_aquint(i,2) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(1) &  a_sims_br(i + sample_half,:) <= quint_temp(2)) > 0);
    alphanz_frac_aquint(i,3) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(2) &  a_sims_br(i + sample_half,:) <= quint_temp(3)) > 0);
    alphanz_frac_aquint(i,4) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(3) &  a_sims_br(i + sample_half,:) <= quint_temp(4)) > 0);
    alphanz_frac_aquint(i,5) =  mean(alpha_star_sims(i + sample_half, a_sims_br(i + sample_half,:) > quint_temp(4)) > 0);
    
    
end

max_alpha = max( max(alpha_star_sims(sample_half+1:end,:)));


%%


disp(' ');
disp(' --------------------------------')
disp([' Table I, Panel A '])
disp(' --------------------------------')
disp( '                                    BR     RE    ControlCost  ');
disp(['  hand-to-mouth agents (%)     :  ' , sprintf('%1.2f   ', 100*[  mean(htm_sims_br(:)), mean(htm_sims_re(:)), mean(htm_sims_tr(:))   ]   )  ]);
disp(['  beta_1 (univariate reg)     :  ' , sprintf('%1.4f   ', [  chtm_reg_br(2), chtm_reg_re(2), chtm_reg_tr(2)] )  ]);
disp(['  beta_2 (multivariate reg)   :  ' , sprintf('%1.4f   ', [  chtm_reg_br2(2), chtm_reg_re2(2), chtm_reg_tr2(2)] )  ]);
disp(['  gamma  (multivariate reg)   :  ' , sprintf('%1.4f   ', [  chtm_reg_br2(3), chtm_reg_re2(3), chtm_reg_tr2(3)] )  ]);
disp(['  delta C  (%) (Ganong-Noel)  :  ' , sprintf('%1.4f   ', [nanmean(cdelta_bigposs_aterc_perc_br(:,1)),nanmean(cdelta_bigposs_aterc_perc_re(:,1)),nanmean(cstardelta_bigposs_aterc_perc_br(:,1))]) ]  );


disp(' ');
disp(' --------------------------------')
disp([' Table I, Panel B '])
disp(' --------------------------------')
disp( '                                             BR       RE      ControlCost   ');
disp([' Avg MPC                                 :  ' , sprintf('%1.4f   ', [  mean_mpc_br, mean_mpc_re, mean_mpc_tr] )  ]);
disp([' Avg MPC | Q5(a)                         :  ' , sprintf('%1.4f   ', [  mean(mpc_br_aquint_sims(:,5)), mean(mpc_re_aquint_sims(:,5)), mean(mpc_tr_aquint_sims(:,5)) ] )  ]);
disp([' MPC | non-HTM                           :  ' , sprintf('%1.4f   ', [  mean(mpc_nonhtm_br), mean(mpc_nonhtm_re), mean(mpc_nonhtm_tr)] )  ]);
disp([' MPC | HTM                               :  ' , sprintf('%1.4f   ', [  mean(mpc_htm_br), mean(mpc_htm_re), mean(mpc_htm_tr)] )  ]);
disp([' Avg MPC | Policy shock as uncertainty   :  ' , sprintf('%1.4f   ', [  mean(mean(mpc_br_infoshock(end,:))), NaN, NaN] )  ]);


%%

htm_ind_temp = htm_sims_br(2:end,:); 

disp(' ');
disp(' --------------------------------')
disp([' Table G.1, Panel (A) '])
disp(' --------------------------------')
disp( '                                                          all agents  HTM   non-HTM  ');
disp(['   Unconditional prob. learning trap                      :' , sprintf('%1.4f   ', [ nanmean(cross_alpha0(:)), nanmean(cross_alpha0(htm_ind_temp == 1)), nanmean(cross_alpha0(htm_ind_temp == 0))])  ]);
disp(['   Trap first period  / life_length                       :' , sprintf('%1.4f   ', nanmean( trap_Tlifelength))  ]);
disp(['   Prob out of trap ever, after first period in trap      :' , sprintf('%1.4f   ', nanmean( trap_prob_outoftrap_ever))  ]);
disp(' --------------------------------')
disp([' Table G.1, Panel (B) '])
disp(' --------------------------------')
disp(['   Mean alpha_t        | shat(ybar) < targ                 :' , sprintf('%1.4f   ', [nanmean( alpha_temp(cross_alpha0(:)==1) ), nanmean( alpha_temp(cross_alpha0(:)==1 & htm_ind_temp(:) == 1) ),nanmean( alpha_temp(cross_alpha0(:)==1 & htm_ind_temp(:) == 0 ) )])  ]);
disp(['   Mean alpha_t        | shat(ybar) > targ                 :' , sprintf('%1.4f   ', [nanmean( alpha_temp(cross_alpha0(:)==0) ), nanmean( alpha_temp(cross_alpha0(:)==0 & htm_ind_temp(:) == 1) ),nanmean( alpha_temp(cross_alpha0(:)==0 & htm_ind_temp(:) == 0 ) )])  ]);
disp(['   Mean   |yt - ybar|  | in-trap                          :' , sprintf('%1.4f   ', [ nanmean( ya_dist_abs(  cross_alpha0  == 1)),nanmean( ya_dist_abs(  cross_alpha0  == 1 & htm_ind_temp == 1)),nanmean( ya_dist_abs(  cross_alpha0  == 1 & htm_ind_temp == 0)) ]) ]); 
disp(['   Mean   |yt - ybar|  | out-of-trap                      :' , sprintf('%1.4f   ', [ nanmean( ya_dist_abs(  cross_alpha0  == 0)),nanmean( ya_dist_abs(  cross_alpha0  == 0 & htm_ind_temp == 1)),nanmean( ya_dist_abs(  cross_alpha0  == 0 & htm_ind_temp == 0)) ]) ]); 
disp(['   Prob (shat_{t+1}(ybar) > targ | shat_t(ybar) < targ)   :' , sprintf('%1.4f   ', [nanmean(cross_alpha0_prime( cross_alpha0_lag(:) == 1) == 0),nanmean(cross_alpha0_prime( cross_alpha0_lag(:) == 1 & htm_ind_lag(:) == 1) == 0),nanmean(cross_alpha0_prime( cross_alpha0_lag(:) == 1 & htm_ind_lag(:) == 0) == 0)   ])  ]);
disp(['   Prob (shat_{t+1}(ybar) < targ | shat_t(ybar) > targ)   :' , sprintf('%1.4f   ', [nanmean(cross_alpha0_prime( cross_alpha0_lag(:) == 0) == 1),nanmean(cross_alpha0_prime( cross_alpha0_lag(:) == 0 & htm_ind_lag(:) == 1) == 1),nanmean(cross_alpha0_prime( cross_alpha0_lag(:) == 0 & htm_ind_lag(:) == 0) == 1)   ])  ]);
disp(['   MPC | trap                                             :' , sprintf('%1.4f   ', [nanmean(mpc_temp( cross_alpha0(:) == 1)), mean([mpc_trap_htm(:,1), mpc_trap_htm(:,2)])]   )]);
disp(['   MPC | out-of-trap                                      :' , sprintf('%1.4f   ', [nanmean(mpc_temp( cross_alpha0(:) == 0)), mean([mpc_notrap_htm(:,1), mpc_notrap_htm(:,2)])]   )]);




